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Volterra Systems With More Than One 
Input Port—Distortion in a 
Frequency Converter 


By S. 0. RICE 
(Manuscript received April 19, 1973) 


Consider a nonlinear system, with memory, which has two input ports 
and one output port. It 1s assumed that the system can be represented by 
a double Volterra series. Two results for such a system are stated in Part I. 
The first is a general expression for the sinusoidal components of the 
output y(t) when the two inputs x(t) and x,(t) are sums of sinusoidal 
terms. The second result is an expression for the power spectrum of y(t) 
when x,(t) is a stationary Gaussian process and x,(t) = P cos pt. Part 
II ts concerned with using results from the theory of Volterra series for 
multi-input systems to calculate the third-order distortion in an tdealized 
frequency converter. 


I. INTRODUCTION 


This paper deals with nonlinear, time-invariant systems with 
memory which (7) have more than one input port, and (7) are driven 
by inputs which are essentially sums of sine waves. 

The paper consists of two parts. Part I is concerned with a system 
which has two inputs, z,(¢) and z,(t), and one output y(t). Two results 
for single-input systems are generalized: (2) an expression is given for 
an arbitrary frequency component of y(t) when 2z.(t) and 2,(t) are 
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finite sums of sine waves, and (iz) an expression is given for the power 
spectrum of y(t) when z.(¢) is a stationary, zero-mean, Gaussian noise, 
and z,(t) is a single sine wave. 

Although most of the discussion in Part I deals with systems having 
two input ports, many of the results can be formally generalized to 
systems with more than two inputs. 

Part II is devoted to an example which shows how results given in 
Ref. 1 for a one-input Volterra system can be used to examine systems 
consisting of a single two-terminal nonlinear element imbedded in a 
linear network containing sources. The transformation from a multi- 
input to a single-input system is based upon Thévenin’s theorem 
(see, for example, Anderson and Leon?). The example treated here 
is a frequency converter using a nonlinear capacitor. Particular atten- 
tion is paid to computing the limiting form of the expression for the 
third-order distortion when the signal and pump amplitudes become 
small. 

The procedure we use in Part II is essentially a systemization of a 
procedure used by Gardiner and Ghobrial* to study the distortion 
performance of a varactor frequency converter. As they point out, 
their treatment differs from the linear time-varying analysis usually 
- employed to study frequency converters. It is appropriate to mention 
here that a promising new general method of computing distortion in 
frequency converters has been developed by R. B. Swerdlow.’ His 
method is based upon the use of Volterra series with time-varying 
kernels. 


Part I. Two Input Ports 


When analysis of the type used to study Volterra systems is applied 
to nonlinear circuits having two input ports and one output port, 
some of the simpler results for one-input circuits can be generalized 
in a straightforward way. Here we state two such generalizations. In 
the first, the two inputs are sums of sine waves. In the second, one 
input is stationary, zero-mean Gaussian noise, and the other input is a 
single sine wave. 

The derivations of the generalizations are not given here because 
they consist of rather straigutforward, although lengthy, applications 
of the procedures used in Ref. 1 to deal with the one-input case. 


II. DOUBLE VOLTERRA SERIES 


Let x.(t) and 2,(¢) be the two inputs, and let the output y(t) be 
given by the double Volterra series 
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ee} ive} 1 foe) foe) foe) foe] 
Wi) = Ev [dur PF dun [doa i 


m=0 n=0 —0 
"Jmjn(U1, °°, Um; V1, °°, Un) 1a Cult — a,) af ty(t — v,), (1) 
r= s= 


where the prime on >~’ means that the term m = n = 0 is omitted. 
The product [J is understood to have the value 1 when the number of 
factors (m or n) is 0, and if n, say, is zero there are no v-integrations. 
The kernel gm,» is a symmetric function of w1, +--+, Um and of v1, +++, Un. 

For the inputs that we shall consider, the (m+ n)-fold Fourier 
transform 


Ganzn( fur; Ny Sums Jars Urs tea) = dui::: dun 


Gmina **> ge “Vn) exp c= J (Uiwu1 Ss a od UnWyn) | (2) 


plays an important role. Here Go;o = 0, w = 27f, and Gn;n 1S a Sym- 
metric function of fui, +--+, fum and of for, ---, fon. 

Much as in Ref. 1, the “harmonic input’? method can be used to 
determine G,,,, from the system equations by setting 


ta(d) = 2 exp (jou), 
- (3) 
(0) = X exp (Jed), 


where the w’s are incommensurable, and solving for the coefficient of 
exp [J (wur e+ ++ um + @o1 +: +++ won)é] in the expansion of y(t). 
This coefficient is equal to Gin;n(fui, +++; fumj for, «++, fon). Note that 
if the system output y(£) is applied to the input of a linear transducer, 
the transducer output can also be expressed as a double Volterra series. 
The transducer output function corresponding to the transducer input 
function Gn; n 1s 


FU (wus pate ta be Wen) Gann fui; peeing 2 fon), 


where F'(jw) is the transfer function of the transducer. 

If, say, n is zero and m > 0, Gmso(fui, fuz, *** » fumj) is equal to 
the coefficient of exp [j(wui + wuz +:---+ wum)t] in the expansion of 
y(t) when x,(t) = 0 and z,(é) is given by (8). 

There is a resemblance between the double Volterra series (1) for 
the two-port inputs z(t), z,(#) (Case A) and the single Volterra series 
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for the special input x(¢) = zu (¢) + x(t) (Case B). For Case B, 
: y (t) = » a ye dui: oe ae durés.(us, ae Uk) 
k 
“Il [u(t — ur) + v(t — u,)]). (4) 


When the product is expanded and the symmetry of é; is used, the 
product can be written as 


m=0 
Setting n = k — mand ums, = v, for s = 1, 2, ---, n carries (4) into a 
form which goes into (1) when @min(wi1, -**; Um; V1, °**, Un) is replaced 
by Gm; n(Ua, met) Um; V1, °°"; Un). 


The results stated below for Case A show a similar resemblance to 
the corresponding Case B stated in Ref. 1. For example, when 2. (¢) 
and 2,(¢) are the sums of sinusoidal terms, the expression for a particu- 
lar component in y(t) for Case A can be obtained from the correspond- 
ing expression for Case B by replacing Gnin in Case B by Gmn and 
inserting a semicolon at the appropriate place in the string of argu- 
ments [as in eq. (9) below]. 


III. SINUSOIDAL INPUTS 


When the inputs z,(t) and z,(t) are sums of sinusoidal waves, an 
expression for any particular component in y(t) can be obtained by 
extending the analysis given in Section VI-B of Ref. 1. There, eqs. 
[1, (188), (139), (140) ] [meaning eqs. (188), (139), and (140) of 
Ref. 1] show that if the input x(t) to a one-input system is given by 


z(t) = > P,, COS wrt, (5) 
r=] 


where the w,’s are incommensurable, then the exp [j(Niwi +::: 
+ N,w,)t], N, 2 0, component of y(t) is 


0 ° L P; Q)Nrt2le 
exp L9(Naes bo + Nato) 22-°- 2, TT writ | 


Gol(fdaitiy (Sidi oes Sdvette (Sadi), (6) 
where Gy = 0, (f.)x denotes the string of k arguments f., fo, °*:, fo, 
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and the subscript n on G has the value 
n= > (N, + 21,). (7) 


Here the notation of Ref. 1 has been changed to bring the statement 
of (6) in line with the notation used in the present paper. 

Methods of computing (6) when the G,’s are constants, i.e., are 
independent of frequency, have been considered by several writers 
(see Kroupa® and Sea and Vacroux®). 

To state the generalization of (6) let 


nN 
Cult = se P,, cos wrt, g(t) = > P, cos w,t, (8) 
r=] r=p+1 


where the w,’s are incommensurable, w, = 2rf,,\ = uw + v, and wand vp 
are positive integers. Then the exp [j7(Niw1 +---+ Nywa)t], N, 2 0, 
component in y(t) is 


; C-) °° r P; 2 Nr+2l, 
exp [y(t tot Mood] EE, | ee aya | 


Gasol (fidatis (—faday (Sedvgetes oo 0s Sudnetter (— Sadia 
(ftv Nepretagy ots Admaem (Adal. (9) 


Here Go;o = 0, and if / or N + l are 0 the corresponding arguments do 
not appear in Gn;n. The values of m and n are 


B 


m = >) (N, + 21;), n= > (N, + 2l,). (10) 


r=1 r=p-+1 


The semicolon in the subscript of Gn;, differs in meaning from the 
semicolon used in [1, (139), (140) ]. The notation (f,), is the same as. 
that in (6) and in [1, (169) ]. The series (9) may either converge or 
diverge, depending on the P’s and G’s. 

Changing the signs of w; and f: in (9) carries (9) into the expression 
for the exp [7(— Niwi + Nowe +--+ Nywa)t] component in y(t), 
etc. [see the discussion below eq. (5) in Ref. 1]. When some of the 
w,’S are commensurable, some of the components in y(t) coalesce and 
can be treated by the method used in []1, (6), (7) ]. 

To examine the case in which 2,(¢) contains a de component, let fi 
and w; tend to 0 in (8) and (9). Then P; is the de component of x, (¢) 
and the exp [j(Nawe +----+ Nywa)t] component of y(t) is the result 
of the coalescence (as f1—0) of (¢) the components exp [j7(Niw1 
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+ Nowe + eae + Nyon)t] for N, = 0, 1, 2, sey, (00 and (12) 
exp [j(— Niwi + Now. +-+-+ Nywr)t] for Ni = 1, 2, ---, ©. When 
(9) and (9) with — f1 in place of f1 are summed over the values of Ni, 
the double sum with respect to J; and Ni can be reduced to a single 
sum by setting k = N;+ 2l, and using the binomial theorem. The 
desired component, namely exp [j(Nowe +---+ Naw )t], in y(é) 
when a(t) = Pi+ P2 cos wot +---+ P, cos w,¢ and 2,(t) is given by 


(8) is found to be 
exp [j(Nawe ++ +++ Noord] 2 Le UG ‘H, lav +h) th | 


‘Ginza (O)ey (fa)ngrtn (—Ffa)tgn ogy Adam (-Adal (1D 


Here n is given by (10), and m = k + (Neo + 2le) +--->+ (N, + 21,) 
when wp = 2 and m = k when p = 1. 

Equation (9) can be generalized to the case of three or more input 
ports in a straightforward way. 


IV. %y(é) GAUSSIAN AND 2,(t) = P cos pt 

The case 2,(é) = P cos pt and z,(t) = I(t), where I(t) is a station- 
ary, zero-mean, Gaussian noise having the two-sided power spectrum 
Wr(f), can be handled in much the same way as was the case x(t) 
= I(t) + P cos pt discussed in Section VII-C of Ref. 1. 

The discrete sinusoidal components in y(é) are given by the ensemble 
average 





(y¥)) = Ps cn exp (jnpt), (12) 
where 
som SP sto 
Smear Sas **%s Su3 fo) = 3 PEM Gr aeset nil fis Sis +004 So 
—foy fry fay ++) Su5 (Snfvdotinis (— Safp)e]. (13) 
Here 2rfp = p, &. = 1forn = 0, s, = — 1 forn < 0, and as in (6), 


(Safp)> denotes a sequence of o arguments, all equal to s,f,. As ex- 
plained in connection with [1, (145) ], Q,.[Wr(f’)] denotes a »-fold 
integration with respect to f;, ---, f, with limits + ©. The integrand 
is Wr(fi)---Wr(f,) times the function [in (13) the function is G] 
of fi, ---, f, represented by all of the terms lying to the right of 
Q.LWr(f') ]. QoL W1(f’) ] denotes the identity operator. 
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The two-sided power spectrum of y(t) is 


Wilf) =X  lenl?8(f — mf) 


Gel 19) >a 6(f — fi —---— fx — nfo) 


°° "(P/gyteeial } 
& lo £ [ny lot Smee iat Rs fol . (14) 


+> 


Replacing Gn.n by Gnyn in eq. (13) for S and substituting in (14) 
gives the expression [1, (175) ] for W,(f) in the single input port case 
x(t) = I(t) + P cos pt. There is a corresponding similarity between 
the one input port formula [1, (16) ] when the two-port expression 
(14) for W,(f) is written out. 


V. EXAMPLE—COMPUTATION OF Gist 


Consider the circuit shown in Fig. 1. The admittance H(f) is linear 
but the resistor R and capacitor C are nonlinear. The voltage across 
R is 
and the capacitance of C depends upon the charge Q, the capacitance 
being a + bQ. The output of interest is the voltage y(¢) across C: 


Q = (a+ bQ)y, 


In + I, = I = dQ/dt. ae) 
The current J,(t) is given by . 
1) = f° nwa — uw) — ye - way (17) 


where h(u) is the Fourier transform of H(f). 
Elimination of J, leads to the circuit equations 


ty = al, + BIZ + y, 
dQ/dt = Ty + / " h(w)Lolt — u) — y(t — udu, (18) 


Q = (a+ bQ)y. 


The G’s corresponding to y can be obtained from (18) by the harmonic 
input method. In using this method it is convenient to work with the 
notation 2, = exp (jwzt) where the w’s are incommensurable. 

In order to get Gio (f1;) we set ty = 21, % = 0, y = ¢1@1 + higher 
harmonies, J, = 732, +---, and Q = qizi +---. The harmonic input 
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Z ae ao with input voltages zr. (t), 2»(t), output voltage y(t), and nonlinear 
and C. 


method states that G1;0 (f1;) 1s equal to c;. Substituting in (18) and 
equating coefficients of 21 gives 


1 = a1 -+- C1, 
jorg, = 11 — A(fi)er, (19) 
di = Ci. 


Solving for c, gives Gio (f1;). Similarly, starting with 7, = 0 and 
Ly = 21 gives Go.1(;f1). The results are 


ties 1 
Gut) = | caries |), 


aa aH 
Go;1 (51) eA eeraneee ie 


where the subscript f; means that f = fi is to be substituted in 
w = Qrf and in H(f). 
To get Gi;i1(fijfe) we set tu = 21, 2» = 22, and assume 


(20) 


Y = Cye1 + Coho + Croke +-°-, 
q, 4121 + 12% + U12%1%2 +-°-, (21) 
Q = 121 + Qo%e + Qiek122 +°°°. 

When (21) is substituted in the circuit equations (18), the coefficients 


of 2, give the equations (19) and therefore c; = G1,0(f1;). Similarly, 
Co = Go.1(;f2). The coefficients of 2:22 give a set of equations which, 


upon solving for ci and using q1 = a¢1, 2 = G2, 42 = — C2/a,t1 = °° +, 
give 
(B/a)(H + jwa);, — J(wi + wee | 
65s = 2c. | ee is 22 
r as (1 + aH + jwaa) 7,44 oe 


Replacing cic2 by G1,0(f1;)Go:1(;f2) gives the required expression for 
G1;1(fij fo) = Cie. 

To get Go.o(f1,f2;) we start with x, = 21 + 22, 2%, = 0 and again make 
the substitutions (21) in the circuit equations (18). The coefficients of 
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2122 give the same set of equations as before because xz, and 2, ap- 
pear only linearly in the circuit equations. We have ¢1 = Gi;o0(/13), 
ty = LA (fi) + goia], gi = aci, and ge = ac, as before, but now 
Co = Gy;o(fe;), 2 = CoLH (fe) + jwoa], and eis = Go;o(fi, fe5). 


To sum up, we have 


Gizi(fi; fe) = 2G1,0( f13)Go;1(; fe) 
X [expression in brackets (22)], (28) 


where Gi,» and G.1 are given by (20). Expressions for Ge.o(f1, f2;) and 
Go:2(;f1, fe) are obtained by replacing the product G1,0Go;1 in (23) by 
Gi,0G1,9 and Go.1Go-1, respectively, and changing the bracket slightly. 
To get Gi,2(fi; fe, fs) we set 2%. = 21, V = 22 + 23 and proceed 
along the lines used to get G1,1, and so on. 
As an example of the use of (23), suppose that 2, = Pi cos wit and 
2, = Pecoswet. Then the exp [j(w: + we)t] component in y is, 


from (9), 
exp Lilor ton)t]] Pais tite] 8 


Similarly, the leading terms in the series for the components of fre- 
quency 2f; and 2f. are given by Go.o(fi, fi;) and Go;2(;fe, fe), 
respectively. 





Part II. Analysis for a Simple Frequency Converter 


Here the circuit shown in Fig. 2 is used as an example to show how 
a multi-input system can sometimes be analyzed by the single-input 
formulas of Ref. 1. The output of interest is the voltage y(t) across the 
nonlinear capacitor C. Thévenin’s theorem is used to replace the circuit 
of Fig. 2 by that of Fig. 3, and a recurrence relation is derived for the 
corresponding G,,’s. The results are used to get an expression for the 
third-order distortion when Fig. 2 is regarded as the circuit for an 
up-converter. 





Fig. 2—Frequency converter with nonlinear capacitor C. 
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VI. REDUCTION OF A MULTIPLE-INPUT SYSTEM TO A SINGLE-INPUT 
SYSTEM 
The system considered here and in the following sections is shown 
in Fig. 2. The admittances H.(f), H.(f), Hi(f) are linear and C is 
the nonlinear capacitor used in Section V. The charge on C is Q(t), the 
current J = dQ/dt flows into C, the capacitance of C is a + bQ, and 
the voltage y(t) across C' is related to Q(t) by 


Q = (a + bQ)y. (25) 


Py is a biasing de voltage. 
The problem is to determine the components of y(t) when 
tu(t) = Py cos wit + Pe cos wel, (26) 
Ly(t) = Py COS wot, 


and w1, we, and w, are incommensurable. 


As far as y(t) is concerned, the analysis of Fig. 2 can be reduced to 
that of the simpler circuit shown in Fig. 3. To accomplish this we 
apply Thévenin’s theorem to the portion of Fig. 2 lying to the left 
of the terminals of C. As far as the exp (jwit) components of y(t) and 
I(t) are concerned, this portion of the system can be replaced by an 
admittance A (fi) = Hu(fi) + Ay(fi) + Ai(fi) in series with the 
(open-circuit) voltage 


; ( H, ) 
— efoit ff —___— ; 
2 A,+H,+ A, fy 


Similar consideration of the remaining components shows that J(t) 
and y(t) can be computed from the circuit of Fig. 3 in which 


H(f) = Hu(f) + Af) + Auf), 


x(t) = poPo + piPi cos (wit + 91) + p2P2 cos (wot + ¢2) 
+ ppPp COS (wet + %), 











= H,(0) Ser cas H(f1) 
po = HO) ’ pier! = Hf)’ (27) 
Hu (fe) fine A, (fo) | 


ig2 = 
tee ay. RE Se) 
The equation for y(¢), namely 


To, = f h(u)[a(t — u) — y(t — u) Jdu, (28) 


where A(u) is the Fourier transform of H(f), can be obtained by 
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he 





| 

| 

y 

! 
Ne 


Fig. 3—Thévenin equivalent of Fig. 2. 


equating two expressions for J, the one on the left being J = dQ/dt 
in which Q = (a + bQ)y = ay/(1 — by). 

It is convenient to subtract out the de components of x(t) and y(t) 
and apply the formulas of Ref. 1 to the portion §(¢) of y(t) which tends 
to 0 when @(t) > 0, €(t) being the time-varying component of «(¢). 
Therefore, in the system equation (28), we make the substitutions 


x(t) = x + &(0), 
y(t) = Yo Si g(t), 


where 2 is the de value of x(t) and yo is the (de) value of y(¢) when 
x(t) = 2. From (27) 20 = poPo, and substitution in (28) gives 
0 = H(O) (#0 — yo). Assuming H(0) + 0 gives yo = Xo = poPo. It 
should be noticed that the substitutions (29) are not strictly necessary 
because the de component in x(t) could be handled (at the cost of 
more work) by the analogue of (11). 

Subtracting the result of substituting x) and yo in (28) from the 
result of substituting (29) in (28) and using 


(29) 








avyot 9) _ ayo _ ag (30) 
1—biy+9) 1—by 1-49 
shows that the system to be analyzed by the single-input formulas of 
Ref. 1 is described by the equations 


eT = f° kee - w — 9 — wd, (31) 


&(t) = piP1 cos (wil + ¢1) + peP2 cos (wot + 2) 
+ ppPp COS (wet + Yp), (82) 


where 
& = a/(1 — by)’, b = b/(1 — byo). (33) 


Here #(¢) and g(t) play the roles that x(t) and y(¢) play in Ref. 1; 
and in the remainder of this paper the G,’s will refer to #(¢) and 9(d). 
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VII. CALCULATION OF THE G,,’S 

The functions Gi(f1), Ge(f1,f2),--- can be computed from (31) by 
the harmonic input method. A guide to the work is furnished by the 
resemblance of our problem to the one described by Fig. 3 of Ref. 1 
and eqs. [1, (42), (48), (106) ]. Expanding the left side of (31) as 


fy abl g() (34) 
l=1 


carries (31) into the form of [1, (106) ] except for the operator d/dt. A 
procedure similar to the one used to deal with [1, (106) ] gives 


H — Qjwd 
an) = (eau), XO = way had 
Go(fi, fo) = bGr(fGr(fo)K (fi + fo), 
Ga(fi, fo, fs) = OGi(fr)Gi(fo)G@i(fa)K (fi + fo + fs) (35) 
‘CK (fi + fe) + K(fi + fs) + K(fo + fs) + 3], 


and the recurrence relation 
Galfiys++, fn) = 3K (fi +--+ + fr) 2 OGD (fay: *+) fn). (86) 


The G’s are the G,’s for the Volterra series for [ y(t) ]’, and formulas 
for computing them are given in [1, (24) to (29) ]. Equation (36) is a 
recurrence relation because G can be expressed as the sum of products 
of Gi, Ge, «++, Gna. By starting with Gi(f1), the relation (86) can be 
used to compute Gs, G3, --- in succession. In the next section, (36) will 
be used to compute G4. 


VIII. COMPONENTS OF y(t) OF FREQUENCY fi + fp AND 2f1 — fe + fo 


In this section we use (6) and the recurrence relation (36) to derive 
expressions for the exp [j(w1 + w,)t] and exp [j(2w1 — we + wp)t] 
components of y(é) in Fig. 2 when (2) Pi, P2, Py are small, (77) fi and 
fe are nearly equal, and (77) H(f) is zero except for frequencies lying 
in narrow bands about the values 


0, fi; fo, fu (37) 


where f, denotes the upper sideband frequency f1 + fp. 

The component of frequency 2f; — fe + fp represents a typical 
third-order distortion product in an up-converter when f, and fe are 
signal frequencies, f, the pump frequency, and fu = fit fp, fo + fo 
the output frequencies. 
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TABLE I—NoTATION FoR VARIOUS VALUES OF H(f) AND K(f) 


Frequency, f H(f) K(f) 
0, fi — fe Ho 0 
Sy fe, 2f1 — fe Ay, Ki 
py Ji 2+ Dp Hy Kp 
fit frol=fu), 2fi — fe t+ fo Hy Ku 
outside bands 0 —2 


As mentioned in Section I, the procedure we use here can be regarded 
as a systemization of a method used by Gardiner and Ghobrial® to 
study the distortion performance of a varactor frequency converter. 
In our notation, the problem they solve is that of determining the 
exp [j(2w1 — we + w,)t] component of the charge Q(t) from the 
system equation 


Vi) = 00 +3Q2 + [ "EG Sana (38) 


Here V (t) is thesum of three sine waves [just as £is in (32) ], J = dQ/dt, 
and k(u) is the Fourier transform of the linear impedance Z(f) in the 
Thévenin equivalent of the converter circuit. It can be shown from 
(38) that the G,’s corresponding to Q(é) can be determined by recur- 
rence from 


Gi(fi) = 1/(@ +JoZ)sy 


—b 
ive 
Sarr oad OF jal | ipmttn 
Now we return to our own problem. From the expression (82) for 
the input @(¢) and the leading terms in the series (6) it follows that the 
exp [j(2w1 — we + w,)t] and exp [j(wi + w,)t] components of y(t) 
are, respectively, 


Gi (fay Fn). (39) 


exp {j[(2w1 — we + wy)t + 201 — g2 + J} 
 (pip2ppP iP oP p/32)[Ga(fi, fr, —f2, fo) +++], (40) 


exp {JL (wi + wp)t + oi + op]} (p1ppP iP p/4) 
*[G2(fi, fo) + (95P3/8)Ga(fi, Joi tay Jn) oh? a (41) 


Only one G, term appears in (41) because we shall assume that 
p1P1/ppP, and p2P2/ppP are small compared to one. 

The function G is given by (35), and the remaining problem is to 
compute G, from the formula obtained by setting n = 4 in the recur- 
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rence relation (36): 

Ga(fi, fo, fa, fa) = 2K (fi + fa + fs + fa) 

BGP (Sa, fay far fd). (42) 


l= 


As explained in Ref. 1 in connection with eqs. [1, (24) to (29) ], we have 


| + 


p G4 (fas fo fay Fx) = (1) (2) (8) (4), (43) 


aa 


p GE (Fay fas fa, fa) = (1) (2) (84) + (1) (8) (24) + (1) (4) (28) 

+ (2)(3)(14) + (2)(4)(13) + (8) (4) (12), (44) 
a GP (fi, fe, fa, fa) = (1) (234) + (2) (184) + (8) (124) + (4) (128) 
+ (12)(34) + (13) (24) + (14) (23), (45) 


oO 


where we have written ‘'(2),’’ for example, for Gi(fe), ‘(34)’ for 
Go(fs, fa), “(234)” for Gs(fe, fs, fa), and so on. 

The next step is to compute the right-hand sides of (43), (44), (45) 
from the expressions (35) for Gi, Gz, Gs; when only frequencies in the 
bands indicated by (37) are allowed to flow. To aid in this, we introduce 
the notation shown in Table I for the values of the function K(f) 
= — 2jw4/(H(f) +jwd) needed for the various G.’s and G;’s. We 
make the usual assumption that the admittance H(f) remains constant 
in each band. 

We consider first the Gi(fi, fi, — fe, fp) in expression (40) for the 
exp [7j(2w1 — we + w,)t] component of y(t). Equation (48) gives 


G2 (fay fas — Fas Fo) = 24G7(fr)Ga(— fo) Gr (Fo), (46) 


where, from (35), Gi(f) = [H/(H + jwd)];. In eq. (44) for Gf 
“(34)” now means G2(— fe, fp) and from the expression (35) for Ge 
and Table I we get 

(34) = 6G1(— fr)Gi(fp)K(— fo + fr) = 6Gi(—f2)Gi (fr) (— 2). 


Similarly, ‘(24)’ means G2(f1,f,) and using the notation K (fit fp) = Ku 
gives 
(24) = bGi(fi)Gi(fy) Ku, 


and so on. Going through all six terms for G{ in this way carries 
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(44) into 
1 é 
31 GP (fi, fis — fo, fo) = bGi(fG1(— fo) Gi (fp) 
-[-2+4+ Kk, +0+ K, +0 — 2). 


Going through all seven terms in G{” carries (45) into 


GP (hs fis — fe, fp) = G3 (f1)G1(— fe)Gi(fr)[Kp(Ku + 1) 


+ (— 2)(— 2) + (0) (Ku) + (Ku) (0) J. 


Substitution of the values of GP (fi, fe, —fe, fr), | = 2, 3, 4, in the 
expression (42) for Gs and combining terms leads to 


Ga(f1, fi, —fey fp) = BGF(S1)Gi(— f2)Gi (fp) 
-K,[2(K, + 1)(Ku + 1) + Ki]. (47) 


When this is put in (40) we get the approximation we have been 
seeking for the third-order distortion (exp [7(2w1 — we + w,)t}) com- 
ponent of y(t). 

The procedure used to obtain (47) can also be used to show that 
the G; in the expression (41) for the exp [7 (w1 + w,)t] component of 
y(t) has the value 


Ga(Fy for for — fo) = BGx(fa)Gi(Fo)G(— fo) 
‘Ku[2(Ki + 1)(Ku + 1) + Ky]. (48) 


If we assume that the two series (40) and (41) converge at about the 
same rate, we can use (48) to get an idea of how large P, can be before 
the leading term in (40) ceases to be a good approximation to the 
typical third-order distortion term. Thus, we expect the leading term 
in (40) to be a good approximation as long as the ratio 


$| PrP pbGi (fy) |*[2(Ki + 1)(Ku + 1) + Kp] (49) 


of the first two terms in (41) is small compared to unity. 

Note that setting f2 = f; and then interchanging f; and f, in the 
expression (47) for Ga(fi, fi, —fe, fp) carries it into the expression 
(48) for Ga(fi, fo, fp) —fp). This is to be expected since Ga(fi, fo, fs, fs) 
is a symmetric function of fi, fo, fs, fa. 
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Picture Coding: The Use of a Viewer 


Model in Source Encoding* 


By J. 0. LIMB 
(Manuscript received March 22, 1973) 


A method is suggested for inserting viewer criteria directly into coding 
algorithms; any complex visual model may be used. The technique is 
applied to a DPC'M-type coder, and a number of variations are compared 
on the basis of entropy, quality, and complexity. It is found that, using 
a simple one-dimensional filter model, the first-order entropy of the DPCM 
signal can be reduced by 30 percent for a high-detail picture with only a 
small reduction in picture quality. Furthermore, by means of a single 
threshold control, one can efficiently trade off bit-rate and picture quality 
over a large range for use in adaptive strategies. 


I. INTRODUCTION 


In early work in picture coding, Graham stressed the role of the 
viewer and Powers and Staras concluded that if large reductions in 
bit-rate are to be achieved they must come from ‘‘nonstatistical”’ 
(perceptual) redundancies.!:? However, there have been few attempts 
_ to explicitly incorporate the viewer in the encoder design. Unfortu- 
nately, there is no general method for handling complex viewer fidelity 
criteria, especially when one is concerned with how pleasing a picture 
appears.t Nevertheless ad hoc techniques have been proposed and 
evaluated and have achieved a certain measure of success.*~® 

Source encoding, in its most general form, can be diagrammed as 
shown in Fig. 1. The first stage is an irreversible operation which 
generates a discrete signal as a result of a quite general multidimen- 
sional quantization process. The resulting discrete signal may still be 
redundant due to the presence of statistical dependencies; these are 
removed in the second stage of reversible processing in which a digital 


a Presented, in part, at the 1972 IEEE International Symposium on Information 
eory. 
T See Ref. 9 for a discussion on viewer fidelity criteria. 
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DETERMINED TO 
CHANNEL 


SIGNAL IRREVERSIBLE REVERSIBLE ENCODER 
SOURCE CODING STAGE CODING STAGE 


Fig. 1—General source-encoding model. 


[perenm 


sequence is assigned to the output of the first stage. Thus, in the 
first stage the properties of the receiver together with the signal sta- 
tistics are incorporated into the quantizing process so that the resulting 
signal just meets the required quality. 

At the output of the first stage, picture quality is established and 
a discrete entropy can be measured. The actual transmission rate will 
then approach the entropy depending on how well the second encoding 
stage is designed to fit the statistics of the source. 


1.1 Receiver-Model Coding 


Algorithm: Components of a picture signal are estimated by some 
method. A test is made to see whether the estimate is adequate by 
testing the estimate on a model of the receiver. If so, the receiver is 
told Gmplicitly or explicitly) that the estimate is adequate. If not, a 
component is transmitted so as to meet the required criterion. 

This type of algorithm will be referred to as ‘‘receiver-model” coding. 
Obviously, it is a rather general approach which can be appended to 
a larger number of existing algorithms; for example, the interpolators 
and predictors summarized by Kortman." In this study we are in- 
terested in applying it to the differential quantizer (DPCM coder) 
although even here it can be applied in many ways. 

In designing a coder to incorporate properties of the human observer 
the most important subjective effect is probably the large decrease in 
visual sensitivity that occurs adjacent to a change in luminance.!— 
An attempt to design a coder based on this effect leads to some form 
of the familiar differential quantizer (DPCM coder).*7:8 

Probably the second most important subjective effect is the change 
in visual sensitivity with average luminance (Weber effect).!4 How- 
ever, in the television situation nonlinearity between applied voltage 
and output luminance in most displays partially offsets this change in 
sensitivity, so that, roughly speaking, noise on an electrical television 
signal is nearly equally visible throughout the luminance range.15.16 

Probably the third most important subjective effect is the spatial 
filtering of small-amplitude, luminance perturbations. It is this third 
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0 VISUAL THRESHOLD YES 
FILTER DECISION NO 
SPATIAL AND 
TEMPORAL 





RANDOM 
NOISE 


Fig. 2—Simple model of visual threshold filtering. 


receiver property which we will attempt to capitalize on through the 
use of receiver-model coding in this paper. 

A simple model that is reasonably successful in explaining the visi- 
bility of liminal stimuli is shown in Fig. 2.131718 Because we are dealing 
with very small perturbations (at least at the neural level) we will ignore 
nonlinearities. The input stimulus on which the model is developed 
is here a small luminance perturbation on a uniform background. 
The stimulus undergoes temporal and spatial filtering and in the process 
is corrupted by noise, represented as an additive random component. 
The filtered signal with the perturbation is compared with the filtered 
background signal. If the difference exceeds a certain threshold then 
the perturbation will be visible.* The model is quite accurate for 
variously shaped stimuli presented on a uniform background with the 
exception that if the stimulus is long (subtended angle >1 degree) 
and thin (subtended angle <5 minutes), it will be significantly more 
visible than the model predicts.!918 . 

The situation is more complex in the case of normal picture evalua- 
tion. First the perturbation is not presented against a uniform back- 
ground and second the perturbation is not directly presented to the 
viewer ; instead it is the difference between the coded picture and the 
viewer’s memory of the original. Thus, although we will use this 
particular filter model it should be upgraded as we understand more 
about the visibility of perturbations in a complex scene. 

In this study we will only be concerned with the spatial effect of 
the visual filter; different shapes have been postulated for the spatial 
impulse response and in one study the Gaussian function was found 
to fit as well as any.13+ However, as we shall see, the performance of 
the algorithm is not sensitive to the exact shape that is used. The 
degree of spread, compared with the size of a picture element, is shown 

* Because of the linearity assumption it does not matter if we filter the difference 


(error) signal or filter the two signals separately and then subtract. 
See also recent work of Ref. 20. 
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Fig. 8—Spatial impulse response of vision: visual point-spread function for 
a Picturephone®-type display reviewed at 36 inches. 


in Fig. 3 for a Picturephone®-type display at a standard viewing 
distance of 36 inches. One should note that this filter is only appropriate 
to threshold vision; once a perturbation is much above threshold it 
may no longer be applicable. 

Note that the efficacy of the filtering operation depends very much 
on viewing distance. Thus, one would expect that at smaller viewing 
distances the eliminated components would no longer be subliminal 
while at larger viewing distances the threshold filtering process could 
be taken further. 


1.2 Coding Algorithms 


Receiver-model coding will be applied to the differential quantizer 
by means of an interpolative algorithm.’ Consider that sample 7 
(Fig. 4) is the last nonzero sample that has been quantized and that 
sample i + j is now being processed. The difference Xi4; — X; is 
formed (where X; is the differentially quantized value of X;), it is 
quantized, and the discrete value of Xi,;, Xi; is evaluated (i.e., 
normal differential quantizer operation). Interpolated values of the 
intermediate samples X41, +++, Xi4j;-1 are then formed from X; and 
X;,; and the error sequence associated with the interpolated values 
is calculated. This error sequence is than processed by the filter- 
threshold circuit to determine whether the errors are visually ac- 
ceptable or not. If the error sequence associated with sample 7 + 7 
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passes the test, the algorithm steps to sample 7 + j + 1 and no new 
value is transmitted. If the test fails, the run is terminated, that is, the 
quantized difference associated with sample z + j — 1 is transmitted. 

There are two distinct forms which the coding algorithm may take; 
free-running or grid. In the free-running algorithm a maximum length- 
of-run is specified in advance for practical reasons. If the interpolation 
attempts to continue beyond the maximum length, a new sample is 
taken and a new run commenced. In most studies the maximum 
length-of-run is 10 pels. In the grid algorithm a fixed set of pels (grid 
elements) is always transmitted and interpolation or extrapolation is 
only applied to the intervening elements. Fixed patterns corresponding 
to every second or every fourth element along a line have been studied 
and the pattern is offset (staggered) from line to line. Grid algorithms 
are studied because in some forms they are very much simpler to 
implement. 

Section II gives the experimental details and describes the basis for 
comparing different algorithms while Section III describes the opera- 
tion and performance of a free-running interpolative algorithm and 
explores the effects of error filtering. In Section IV we describe and 
compare the operation of a number of different grid algorithms in- 
cluding one which involves but a minor modification to the normal 
differential quantizer. 


II. EXPERIMENTAL ARRANGEMENT 


The different algorithms were evaluated using a computer facility. 
The 8-bit digitized pictures are read from a digital disk, line at a time, 






AMPLITUDE =——> 


i i+ i+j-1 i+j 
DISTANCE IN PICTURE ELEMENTS 


Fig. 4—In description of an extrapolative threshold coder. 
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processed, and then stored in a digital frame store for direct viewing 
on a television monitor. The picture consists of 250 lines with 210 
elements in each line. The picture is generated and displayed as a 2:1 
interlaced picture at 30 frames (60 fields) per second; hence adjacent 
lines in the picture originate in different fields. This format is similar 
to the Picturephone format. 

In evaluating the picture we look at a single frame, repeated at 
30 frames per second; thus temporal effects are not considered. The 
picture quality is slightly better when viewing a “‘frozen’’ frame of a 
differentially quantized picture since ‘‘edge busyness’”’ and certain 
random noise components are noticeably less objectionable in the 
frozen situation contrary to the findings for white noise.”! 


2.1 Differential Quantizer 


The normal differential quantizer is the vehicle with which the 
various algorithms will be tested. The 13-level companded quantizing 
characteristic is given in Table I. The differential quantizer has no 
integrator “leak’’ but the integrator is reset at the beginning of each 
line. 

The results will be given mainly in terms of two different pictures. 
The first picture is the familiar ‘‘Karen’”’ which by most measures 
would be regarded as active and is fairly difficult to code if both the 
soft hair and the sharp stripes are to be preserved. The second picture 
is much simpler having a large flat background and is referred to as 


TABLE I—QUANTIZER CHARACTERISTIC OF 13-LEVEL 
DIFFERENTIAL QUANTIZER 
(expressed in 1/128ths of the p — p amplitude) 


Level Number Decision Level Representative Level 

0 0 
1 

+1 2 
3 

+2 4 
6 

+3 8 
11 

+4 14 
18 

+5 22 
27 

+6 32 
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“Lamp.” A third picture (‘Birdcage’) is occasionally used; it is 
intermediate in complexity between the previous two. 

The picture quality of the differentially quantized signal is only 
distinguishable from the 8-bit digital signal by careful comparison; 
there is a slight increase in background noise and very small amounts 
of slope overload and edge-busyness. The discrete, first-order entropies 
of the three pictures after coding by the differential quantizer are 
3.10, 2.79, and 2.37 bits/pel for Karen, Birdcage, and Lamp, respec- 
tively. The second-order entropies are 2.92, 2.61, and 2.20 bits/pel, 
respectively. Thus, little would be gained in the second stage of coding, 
the reversible stage, by any attempt to remove higher-order 
redundancy. 


2.2 Quality 


One difficulty in documenting the performance of coders lies in 
specifying the quality of the processed pictures. 

One can divide picture quality into different ranges by using a set 
of criteria. Consider the following three: 


1. Difference just detectable by a skilled observer between the 
processed and unprocessed pictures in an A~B comparison with 
no restriction on viewing distance. 

2. Defects just noticeable to a skilled observer at standard viewing 
distance (36 inches approximately 7H) for a picture with which 
the observer is familiar. 

3.* Defects just noticeable to a skilled observer, at standard view- 
ing distance when the observer has no knowledge of the original 
picture. 


The picture quality of criterion 1 is probably the most frequently 
used ad hoc criterion but it is unnecessarily severe for visual communi- 
cation purposes and, if employed, would result in a significant in- 
crease in bit-rate over that required by criteria 2 and 3. In this study 
the author has attempted to specify the qualities of coded pictures 
using criteria 2 and 3. This is inevitably an approximate process and 
as a consequence a range is given rather then a specific value. Ap- 
proximate as this process is, if it enables a rank ordering of coding 
strategies it will have served its purpose. 


* Where the viewer was familiar with the test picture a conscious effort was made 
to disregard defects that depend on knowledge of the original picture. For example, 
noticing a loss of fine detail in the hair region of ‘‘Karen’’ depends on memory of the 
original; noticeable slope overload on the other hand generally appears as an un- 
natural distortion. 
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2.3 Bit-Rate Calculation 


Picture quality and bit-rate are the two vital measures of coder 
performance. In this study we are concerned primarily with the first 
coding stage of Fig. 1, the multidimensional quantizing stage. But 
the final bit-rate will also depend on how thoroughly the second stage 
is implemented. However, what we will do is to calculate entropies 
of the signal after the first stage of coding, the rationale being that the 
figure represents a bound on what is obtainable in practice. In some 
instances variable wordlength coding, with buffering, will yield a 
data rate that is within a few percent of the entropy figure.??:23 In other 
instances more complex coding will be required to approach the entropy 
figures, particularly for source alphabets which contain a highly 
probable event where something akin to runlength coding would be 
required. 

The performance of the algorithms has been assessed by calculating 
the entropy under the assumption of two different types of reversible 
encoding. They are:* 


Code I. All pels in the run are processed in the same way (with the 
same code). This is the simplest but most inefficient method. 
The bit-rate bound is obtained by calculating the first-order 
entropy of the signal; 


N 
H, = — 2 pi log pi, 


where p; is the probability of occurrence of each event (a 
quantizer level or an interpolate command) and N is the 
total number of different types of event. 

Code II. A separate code is used for each run position. That is, the 
first element in a run uses code 1, the second element in the 
run uses code 2, etc. The run is terminated by the sampled 
pel. The entropy is then given by 


M 
Hy = 2 hq, 
j= 


where h; is the entropy of events in the jth position of the 
run and q; is the probability that an event will be in the jth 
position of a run. 


*In some of the algorithms to be discussed the information indicating that an 
element has been successfully estimated at the transmitter is obtained indirectly 
by the receiver from the coded bit stream. In other algorithms an additional code 
word is appended for this purpose. 
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Fig. 5—Variation of entropy with the position of the element in the run; free- 
running interpolative algorithm with a maximum runlength of 10. Subject—Karen. 


The entropy of the signal changes significantly depending on the 
position in the run. (This is shown in Fig. 5 where the first-order 
entropy of the differentially quantized signal is plotted as a function 
of the position in the run for a free-running interpolative algorithm 
having a maximum runlength of 10.) This change in entropy is ex- 
ploited in code II (but not code I). Where the average length of a 
run is large, a practical realization of a code II coder could well result 
in a type of runlength encoding. 

In summary, Hi, can be regarded as the lower bound on data rate 
when each element is coded in the same way while H: is a lower bound 
when run contiguity is exploited. 

There are a great number of different techniques for reversibly 
coding the discrete output of the first coding stage (Fig. 1); by speci- 
fying the abovementioned two entropies we can concentrate more on 
the irreversible stage without getting overly involved in exactly how 
the second-stage coding will be achieved. The entropies are always 
given as bits/active (or unblanked) picture element. 


III. RESULTS : FREE-~RUNNING ALGORITHM 


The details of the interpolative algorithm are summarized in the 
flow diagram of Fig. 6. Bookkeeping operations like entering a new 
line, testing for the end of a line, and gathering statistics are not 
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Fig. 6—Flow diagram for the element processing of the free-running interpolative 
algorithm. J denotes the last element in the previous run, J denotes the current length 
of the run being processed, and J + J denotes the element being currently processed. 


shown. We will first discuss (Section 3.1) the efficiencies obtained with 
the two methods of reversibly coding the discrete output. Neither the 
shape of the filter function nor the maximum length which the algo- 
rithms can run before a new run is forcibly commenced is varied in 
the above comparisons. The effect of varying these two parameters is 
described in Sections 3.2 and 3.3. Some observations are made on 
free-running algorithms in Section 3.4. 
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Fig. 7—Entropy of the free-running algorithm as a function of the threshold 
The measurements of entropy are made under the assumption of two different types 
of reversible code. The performance for the codes is very similar. 


3.1 Comparison of Reversible Coding Methods 


Figure 7 summarizes the results obtained by applying receiver- 
model coding interpolatively to the 13-level differential quantizer. * 
For computational simplicity, the filter used in this case has a rec- 
tangular impulse response three elements wide (i.e., corresponding to 
an average over three elements). 

As the threshold is raised on the filtered error sequence, more and 
more elements are interpolated. Consequently probability distribu- 
tions become more peaked and the entropy drops. At the same time 


“The “relative” threshold is, in fact, one-fifth the threshold value, in 128ths, 
applied to the filtered error signal. 
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Fig. 8—(a) Karen—processed by normal 13-level differential quantizer, 1st order 
entropy 3.10 bits/pel. (b) Picture processed by free-running algorithm, 2.0 bits/pel; 
picture quality is criterion 3 or worse. (ec) Unprocessed picture of ‘‘Lamp.”’ 
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Fig. 8 (continued). 


the picture quality is reduced in low-detail areas of the picture as 
soft detail and texture become blurred. Edges and high-detail areas, 
however, remain unaffected until very large thresholds are reached. 

Consider the results for Karen. As the threshold is increased, the 
entropy drops from 3.1 bits/pel with a normal differential quantizer* 
to about 2 bits/pel at which point there is quite noticeable smearing 
in low-detail areas. Also shown on the curves are the criterion 2 and 
criterion 3 ranges (Section 2.2). Not until the threshold is raised to a 
value of 0.9 and the entropy has fallen to 2.4 bits/pel does the change 
in picture quality become visible when compared with a normal dif- 
ferential quantizer, other than by close A-B comparison. The normal 
differentially quantized picture is shown in Fig. 8a while the picture 
coded with a threshold of 1.5 (2.0 bits/pel) is shown in Fig. 8b. 

The results obtained with the simpler picture ‘“‘Lamp” (Fig. 8c) 
are similar to those obtained for Karen except that the advantage is 
somewhat greater; the rate is halved in going from the normal dif- 


*It is necessary to send additional information to explicitly inform the receiver 
when to interpolate and when not to. It is this additional information which prevents 


the entropy of the coded signal from converging to the value of the normal differ- 
ential quantizer. 
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ferential quantizer to the end of the criterion 3 range. It is to be ex- 
pected (see Section V) that low-detail pictures will be more amenable 
to receiver-model coding given the present model. 

There is surprisingly little difference in efficiency between the two 
reversible codes, particularly for Karen where the statistics for the 
highly detailed parts swamp the peaked distributions obtained in the 
low-detailed parts. In such instances an adaptive strategy would be 
of some help.* The complexity associated with implementing the simple 
code (code I) does not change with the maximum permitted length 
of run; for the variable code (code II) there is a proportional relation- 
ship since a code dictionary would need to be stored for each run 
position. Consequently, it is important to know how the entropy 
changes with the maximum length of run that is permitted. For the 
moment we may conclude that unless the more complex codes can be 
implemented simply or that channel capacity is at a premium then 
the simple code is probably adequate. 


3.2 Visual Filter Function 


The psychological literature is replete with different estimates of 
what the shape of the visual point-spread function should be. It was 
hoped that we could add something to the debate by investigating 
different functions in the coding model to see which shape gives the 
best results. In one experiment the shape of the function was varied 
keeping the spread of the function constant; the spread was measured 
by the first moment of the absolute value of the spread function. In a 
second experiment the spread of the filter was varied keeping the shape 
constant. Bear in mind that because our algorithm works only along 
the scan-line we cannot take full advantage of the two-dimensional, 
spatial, point-spread function. Consequently, we should really think 
of a line-spread function, the rationale being that in the worst-case 
situation the stimulus being filtered would have large vertical extent 
and hence the line-spread function would be appropriate. 


3.2.1 Effect of Shape 


Varying the shape of the filter function has little effect on coding 
efficiency (Fig. 9). The filter shape was varied from rectangular to 
quite peaked keeping both the area under the function and the first 
moment of the absolute value of the function constant. The threshold 
is also constant at 1.0. The square, crosses, and dots of Fig. 9 denote 
functions with widths of 3, 5, and 7 pels respectively. The values of 
the functions are given in Table II. For the interpolative algorithm 
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g. 9—Effect on entropy (code II) of varying the shape of the filter function. 
Tm mean of the impulse response is: []—3 elements, X—5 elements, O—7 ele- 
ments. Threshold decisions are very insensitive to the shape of the filter function. 


with a maximum runlength of 10 pels there is an increase in bit-rate 
from 2.32 bits/pel for the rectangular function to 2.41 for the most 
peaked function; any accompanying change in picture quality was too 
small to notice. 


3.2.2 Effect of Spread 


The spread of the filter function, on the other hand, has far more 
effect on the picture quality and entropy than does the shape, as can 
be seen from Fig. 10a. A rectangular function was used and the spread 
was varied keeping the area under the impulse response constant and 
the threshold fixed at 1.0. The picture quality changed from almost 


TABLE IJI—WEIGHTING COEFFICIENTS OF TRANSVERSAL FILTER 
(The filter shape is symmetrical with A being the central element) 


Filter Number A B C D 
1 0.333 0.333 0 0 
2 0.4 0.275 0.025 0 
3 0.45 0.231 0.044 0 
4a, 0.5 0.188 0.062 0 
4b 0.5 0.156 0.062 0.031 
5 0.55 0.103 0.081 0.041 
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SPREAD OF RECTANGULAR FILTER FUNCTION IN PELS 
Fig. 10a—The effect on entropy (code II) of varying the width of the filter func- 


Hon: The overall spread of the function has a strong effect on entropy. Subject— 
aren. 


-- CONSTANT QUALITY 
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0 0.5 1.0 1.5 2.0 
RELATIVE THRESHOLD 


Fig. 10b—Curves of entropy versus threshold for filter functions having spreads 
of 1, 3, 5, and 7 elements for the free-running interpolative algorithm (Karen). 
The dashed curve passes through each of the full curves at points of approximately 
constant picture quality. A spread of between 3 and 5 elements gives the lowest 
bit-rate for standard viewing distance. 


criterion 1 quality with a spread of 1 pel to worse than criterion 3 
quality when the spread was 7 pels. 

An attempt was made to determine the most suitable filter spread 
for a picture having criterion 2 quality (standard viewing distance). 
Figure 10b gives curves of entropy versus threshold for rectangular 
filter functions of different spread. The dashed curve is a line of ap- 
proximately constant picture quality. It was determined by making 
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pair-wise comparisons between a reference picture obtained using a 
threshold of 1.0 and a filter spread of three and pictures from the other 
spread curves. With the filter fixed at a particular value of spread the 
threshold was varied until the picture quality matched that of the 
reference picture. From the figure it can be seen that a spread of 
between 3 and 5 pels gives the lowest bit-rate for the standard viewing 
distance. 


3.38 Effect of Maximum Runlength 


The effect of changing the maximum permitted runlength is shown 
in Fig. 11. Interestingly, there is very little increase in bit-rate as the 
maximum runlength is reduced to as little as 4 pels, particularly for 
code II. Even for the low-detail picture (Lamp) where the average 
length of a run is much longer, the increase in entropy is still small. 
Bearing in mind that code II becomes much simpler to implement for 
short maximum runlengths there appears to be little reason to use long 
runlengths. 


3.5 


3.0 


2.5 


2.0 


RATE IN BITS PER PEL 





MAXIMUM RUNLENGTH OF ALGORITHM 


Fig. 11—Entropy as a function of the maximum runlength for Code I (dashed) 
and Code II (full). Note there is little increase in entropy for Code II as the maximum 
runlength is reduced from 10 to 4. 
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Fig. 12—Pictures showing the effect of changing the size of the picture element 
with the filter function, as measured at the eye, maintained constant: (a) original 
8-bit signal, (b) processed, with threshold = 1.5 and H = 1.38 bits/pel, (c) original 
8-bit signal, 4 lineal size, (d) processed with same threshold as in (b), H = 1.17 
bits/pel. It is the quality difference between pictures of the same size that should 
be compared, not the relative quality of the two processed pictures. 
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Fig. 12 (continued). 


3.4 Discussion 


The preceding experiments suggest two ways for decreasing bit- 
rate at the cost of decreased picture quality. First, it can be decreased 
by increasing the threshold as shown by Fig. 7. Second, it can be 
decreased by increasing the spread of the filter function as shown by 
Fig. 10a. The picture, Karen, was coded to have an entropy (Code IT) 
of 1.81 bits/pel by reducing the quality (lower quality than criterion 3) 
in the two ways described above. For the first method the filter was 
rectangular with a spread of three elements while for the second method 
the filter was again rectangular but with a spread of seven elements. 
Both methods gave similar picture quality with the narrow-filter/ 
high-threshold combination of the first method being, perhaps, slightly 
better. The improvement in sharpness of the first method was partly 
offset by the reduction in granularity and blotchyness of the second 
method. 

If a particular filter, at normal viewing distance, produces a picture 
that is just distinguishable from a high-quality original then doubling 
the spread of the filter function should produce a picture at twice the 
viewing distance which is again just distinguishable from the original. 

I have tried to demonstrate this prediction with Fig. 12 by repro- 
ducing a comparison pair of pictures at half-size to correspond to the 
situation where the viewing distance is doubled. It is the dzfference in 
quality between pairs of pictures at the same viewing distance that 
should be compared, not the comparative quality of the processed 
pictures. 

One factor that could upset such a comparison is that the smaller 
picture has a greater scanning line density. The filter function operates 
in one dimension only and to the extent that deleted picture com- 
ponents are uncorrelated from line to line, vertical filtering taking 
place in the eye will tend to favor the smaller picture. An intuitive feel 
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Fig. 13—Picture of the filtered error signal for the processed picture of Fig. 12b. 


for the correlated nature of the error signal is obtained from Fig. 13, 
in which a certain amount of picture structure is evident. 


IV. RECEIVER—-MODEL CODING WITH GRID ALGORITHMS 
4.1 Introduction 


One can take advantage of the filtering action of vision without 
explicitly filtering the error signal. To appreciate this, let us consider 
the following grid algorithm. Every grid element (sampled point) is 
reproduced with full accuracy (e.g., 7 or 8 bits). The intermediate 
elements (referred to as “conditional points’) are reproduced as the 
average of the adjacent pels, Xis1 = (X; + Xi42)/2, if the error 
(Xis1 — Xsi1) is small (see Fig. 14). Otherwise, the error quantity is 
quantized and transmitted. In determining whether X;,1 is an adequate 
representation of X,41, the error signal adjacent to pel (¢ + 1) must be 
filtered. However, the error at pels 7 and (7 + 2) is virtually zero so that 
for a filter that consists of a three-point average it is only necessary to 
examine the error introduced at pel (¢ + 1). 

Kretzmer®! proposed a coding scheme similar to the above in which 
every fourth pel is always coded with 7-bit accuracy (i.e., 4:1 grid 
algorithm). The intermediate points are estimated by linear inter- 
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polation and the difference between the input and the estimate is 
quantized and transmitted. The midpoint in each quad is quantized 
more accurately than the quarter and three-quarter points. Fukushima 
and Ando*> experimented with a very similar scheme in which every 
fourth point was transmitted with 6-bit accuracy and the intermediate 
points were transmitted using three levels. A final bit-rate of 2.7 
bits/pel was achieved. They also investigated two-dimensional 4:1 
algorithms. Connor has investigated a 2:1 grid algorithm (column 
coder) which uses two-dimensional prediction for differentially coding 
the grid points.2® Pease?’ has applied what amounts to a 2:1 grid 
algorithm between fields of a television picture. All points in one field 
are estimated as the average of the four surrounding points coming 
from the previous and next fields. Only when this prediction breaks 
down is additional information sent about the interpolated field. In 
the presence of movement the four-way interpolation is less accurate 
and the number of pels that require correction increases somewhat. 
Notice that all the above schemes transmit two or more different types 
of amplitude information; the grid points are transmitted absolutely 
(or differentially, relative to one another) while the conditional 
points are transmitted as a correction to the estimation. These schemes 
will therefore be referred to as error transmission schemes. 

In this section we will examine a number of grid coding schemes. 
For the most part they differ from the above schemes in that only one 
type of amplitude signal is transmitted so that all amplitude informa- 
tion is decoded in the same way (direct transmission). The distinction 
is best appreciated by considering a specific example. Take the inter- 


AMPLITUDE 





i i+1 i+2 
DISTANCE IN PELS 


Fig. 14—Definition of locations and values of elements used in discussion of grid 
algorithms. 
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polative algorithm: if pel 7 has already been encoded (Fig. 14), pel 
(¢ + 2) is then encoded differentially from pel 7. From the encoded 
values of pel ¢ and pel (¢ + 2), (X;:, Xi+2), X41 is formed. The error 
signal (Xis1 — Xis1) is tested against the threshold; if it exceeds 
threshold, pel (7 + 1) is differentially coded from pel 7 and pel (7 + 2) 
is differentially recoded from pel (7 + 1). Thus it can be seen that the 
interpolated value X;.1 is only retained when the interpolation is 
adequate; otherwise it is discarded. Furthermore, the quantizing scales 
for pels 7 and (¢ ++ 1) can be the same as for normal differential 
quantization since in high-detail areas the interpolation generally fails 
and each element is predicted from the previous element. In practice, 
a check is made to determine whether slope overload will occur in cod- 
ing pel (4 + 2); if this can happen pel (z + 1) is then coded and pel 
(¢ + 2) is recoded, differentially, from pel (¢ + 1). Thus, in high-detail 
parts of the picture, pel (¢ + 1) is rarely interpolated and the coding 
operation differs little from normal differential quantization. In low- 
detail parts of the picture, where the interpolation process is usually 
adequate, again the coding process is normal differential quantization, 
but with twice the normal sample spacing.* 

Errors will occur at pels 7 and (7 + 2) because differential quantiza- 
tion has been used and these errors will, because of the visual filtering 
action, affect the visibility or the error occurring at pel (7 + 1). 
Hence the encoding will be more efficient if filtering is used. But, as 
we will see, a three-point filter does not differ much from a single- 
point filter because the errors made at pels 7 and (7 + 2) are limited by 
the number and spacing of the quantizer levels and cannot be sub- 
jectively large if adequate quality is to be obtained. 

In comparing the error transmission and direct transmission schemes, 
it can be seen that the decision on whether or not to transmit the con- 
ditional elements is the same in both cases. The error transmission 
scheme has the advantage that the estimate is a better prediction 
than the previous sample, and hence the correction signal, where it is 
necessary to transmit it, will be smaller. However, the disadvantage 
is that since the grid points are transmitted as differences from a 
point two pels away, the amplitude of the differences and hence the 
entropy associated with them will be larger. In practice this will 
increase complexity since the quantizer will need to have more levels 
to handle the larger changes. In Section 4.2 an error transmission 
scheme will be compared with a number of direct transmission al- 
gorithms and it will be seen that there is very little difference in per- 
formance between the two types of schemes. One would expect the 


RECEIVER-MODEL CODING 1293 


performance to converge for low-detail pictures since the number of 
points which are not successfully interpolated becomes very small 
and the encoding of the remaining points is then very similar. 

In the free-running algorithms a special code word was used to 
inform the receiver when to interpolate. For the grid algorithms an 
interpolate command has been inserted in a special manner. On the 
conditional samples only, the zero differential quantizer level is used 
to denote the interpolate command: this means that when the signal 
is not being interpolated the zero level cannot be used; instead the 
signal is forced to take on the next closest level, either the positive or 
negative inner level. This affects picture quality very little since, 
firstly, a zero level is rarely used on the conditional samples and, 
secondly, since interpolation generally fails in the vicinity of large 
luminance changes, the small error introduced by deleting the zero 
level is largely masked by the consequent luminance change. 

Implementation of the grid algorithm becomes even simpler when 
we consider two variations, a modified form of the interpolative (MI) 
algorithm and an extrapolative algorithm. The MI algorithm is quite 
similar to the interpolative algorithm; the next grid point is not 
quantized prior to interpolation. This means that it is only necessary 
to quantize each element sequentially just as one does in normal 
differential quantization (when a pel is adequately interpolated, the 
classifier output is simply forced to a zero prior to processing by the 
local [and distant] decoder and the next element [pel 7+ 2] is 
processed in the normal manner [see Fig. 14]). In the extrapolative 
algorithm the method used to estimate the conditional sample is the 
same as the method of extrapolation for the coding process (Le., 
previous sample prediction) and hence the need for an extrapolate 
command is obviated. The algorithm is then only slightly different 
from normal quantization, especially if the error occurring at the 
conditional sample is taken as the filtered value (the scheme described 
in Ref. 4 under the name “‘Level Variable Sampling Scheme’’). 


4.2 Comparison of Free-Running and Grid Algorithms 


The performance of both a 2:1 and a 4:1 MI, grid algorithm are 
compared with the free-running extrapolative algorithm in Fig. 15. 

The maximum reduction that can be obtained with the 2:1 algorithm 
is a halving of the bit rate. Long before this point is reached the curve 
starts to flatten out and unless very large thresholds are used the 
picture quality remains high. Within the obtainable range of picture 
quality the 2:1 algorithm performs almost as well as the free-running 
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Fig. 15—Comparison of the performance of free-running and grid algorithms. 
The 4:1 grid algorithm performs equally as well as the free-running algorithm. 


algorithm. By going to the 4:1 algorithm, a larger picture quality 
range can be accommodated without going to very large thresholds. 
In the criterion 2 range the grid algorithm seems slightly better than 
the extrapolative free-running algorithm while in the criterion 3 
range the free-running algorithm is marginally better. 


4.3 Comparison of Three Grid Algorithns—Error-transmission, MI, 
and Extrapolative 


Since the MI and error-transmission algorithms are the most alike, 
we will compare them first. The error-transmission algorithm uses a 
19-level differential quantizer. This is obtained from the 13-level 
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quantizer by adding additional outer levels. The filtered error signal 
is obtained by summing the error at the estimation point and the 
two adjacent grid points. The MI algorithm uses the usual 13-level 
quantizer and the filtered error signal is the sum of only two error 
terms. The quantizing error occurring at the grid point to the right 
of the point being interpolated cannot be included since this point is 
not quantized until after a decision has been made on the conditional 
point. 

The white markers in Fig. 16 indicate those conditional points in 
the two algorithms for which the filtered error signal is above thresh- 
old. Hence these points are not adequently represented by the estimate 
(the relative threshold is set at 1.5 for both algorithms). The distribu- 
tion of markers is quite similar, especially when one bears in mind 
that the error summing procedure is different in the two cases. The 
picture quality and bit-rate is also very similar (see Fig. 17), which 
stands to reason since the signal is processed identically in those 
parts of the picture where there are no markers. The algorithms were 
evaluated on other pictures. In each case picture quality and bit-rate 
were very close. 

The extrapolative (like the MI) algorithm uses a 13-level quantizer 
and sums the error over only two pels. The estimation procedure 
(zero-order-hold) is not as effective as linear interpolation and, as a 
result, the number of conditional points that need to be transmitted 
is very much larger for a specific threshold. A consequence is that the 
curve of entropy versus threshold lies above the other curves except 
at higher thresholds. Here, the curves converge since the only condi- 
tional points still being transmitted are edge points. The picture 
quality is not quite as high as that obtained with the other two al- 
gorithms with the defect appearing as a granularity in flat, dark 
regions of the picture. Although the granularity is also present for the 
other two algorithms it is significantly attenuated by the interpolative 
averaging. 


4.4 Effect of Filtering 


As indicated previously, the effect of filtering for the 2:1 grid al- 
gorithm will not be very strong since when the error is evaluated at 
each conditional point the error permitted at the adjacent points, 
which are quantized with full accuracy, will be quite small. Even so, 
there is a small increase in the number of conditional flags that are 
transmitted in going from the single-point filtering to the two-point 
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Fig. 16—Markers showing conditional points that were updated with a threshold 
of 1.5: (a) error transmission algorithm, (b) MI algorithm. 
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Fig. 17—Relative performance of three different 2:1 grid algorithms. The extrap- 
oe algorithm is slightly inferior to the MI and dual-mode algorithms. Subject— 
aren. 


filtering (error at conditional point plus the error at previous point). 
This, in turn, results in a small increase in entropy (from 2.17 to 2.20 
bits/pel). 

For the 4:1 fixed-point algorithm the difference between single- 
point and three-point filtering is larger. The conditional points that 
are transmitted have been marked in Fig. 18 where, for single-point 
filtering, the threshold is 0.9 and the entropy is 2.08 bits/pel and for 
three-point filtering the threshold is 1.5 and the entropy is 2.04 bits/ 
pel. In this case, however, the effect on picture quality is more notice- 
able. With the broader filter low-detail areas are reproduced better 
while medium-detail areas appear more noisy. At normal viewing dis- 
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tances the broad filter is preferable while for close scrutiny the single- 
point filter is better. 

There is no reason why filtering could not take place in two or three 
dimensions in which case more elements would be involved and the 
accuracy with which the picture was encoded could more accurately 
match perceptual requirements for a given viewing situation. 


V. DISCUSSION 


As we have seen, the receiver-model coding algorithm with the simple 
threshold model of Fig. 2 tends to work best on low-detailed pictures. 
There are two reasons for this: (z) In detailed parts of the picture the 
estimation procedure is not as good as in low-detail areas; (77) The 
threshold model, as described, is a simple, low-pass filter model and 
does not incorporate the effects of masking by adjacent signal com- 
ponents such as occurs when an element lies close to a large change 
(spatially or temporally) in luminance.* 

The receiver-model coding concept, as stated, does not depend on 
any specific receiver model. As better models of the human viewer 
are obtained they can be incorporated directly into the encoding 
operation. In essence it is a three-step operation: estimation, testing, 
and, if necessary, more accurate recoding. There is an intrinsic separa- 
tion between the source-property operation (estimation) and the 
recelver-property operation (testing) and as such the technique will 
be suboptimum. Performance could undoubtedly be improved by 
cycling through the estimate-test-recode sequence iteratively.?® The 
interesting, practical question would be, is the improved performance 
worth the added complexity? 

In all the coders described here the bit-rate—picture-quality operat- 
ing point is determined by means of a single threshold control. This 
means that it is a relatively simple matter to dynamically alter the 
operating point in response to some system requirement. An example 
occurs in frame-to-frame coding where the moving area is transmitted 
as an element-differentially-quantized signal. As the buffer fills in 
response to increased movement the threshold is raised so as to keep 
the data-generation rate more uniform.°? 


VI. SUMMARY AND CONCLUSIONS 


Receiver-model coding is a powerful, though not optimal, technique 
for incorporating properties of the human observer into the picture 


*Some practical coding strategies have been developed that take advantage of 
spatial masking effects.+'6 
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encoding process. In essence, components of the signal are estimated 
according to some algorithm. The difference between the actual 
signal and the estimate is processed in a model of the receiver to 
determine if the estimate is adequate. If so, the receiver is informed of 
this; if not, additional information is transmitted to improve the 
estimate. 

The receiver-model coding concept may be applied in many dif- 
ferent ways and the visual model may range from very simple to very 
complex. In this paper I have used the differential quantizer (DPCM 
coder) as the basic vehicle with which to investigate receiver-model 
coding, and the visual model is a one-dimensional low-pass filter. 
Three types of estimation are investigated: extrapolation, interpola- 
tion, and a simplified form of interpolation referred to as ‘‘modified 
interpolation.” It is important to bear in mind that the estimation is 
used to help determine which components need to be transmitted and 
does not indicate how the components are transmitted. In nearly all 
examples considered here the transmitted component is a simple 
difference signal which is decoded by adding the difference to the last 
decoded value. 

Coders are divided into two separate classes, free-running algorithms 
and grid algorithms. In the free-running algorithms the estimation 
procedure may continue in a single run until the estimate fails with 
the proviso that the length of the run may not exceed a specified 
maximum. With the grid algorithm a fixed set of elements is always 
transmitted (e.g., every second or every fourth element). The interest 
in grid algorithms stems from the fact that they are more easily 
implemented. 

The free-running interpolative algorithm gives a reduction in entropy 
of approximately 30 percent for high-detail pictures and 50 percent 
for low-detail pictures for a small loss in picture quality when the 
picture is evaluated by observing a single “frozen” frame on a high- 
quality CRT display. 

Two reversible coding strategies were explored for converting the 
quantizer output to a binary code. Code II gives an advantage of 
between 0.1 and 0.15 bits/pel over Code I when using a maximum 
runlength of ten elements; the relative advantage of Code II over Code 
I about doubles when the maximum runlength is reduced to four 
elements. 

The effect of the threshold filter function on the coding operation 
was explored by varying the shape of the filter function while keeping 
the spread of the function constant and then, in a second experiment, 
keeping the shape constant and varying the amount of spread. While 
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the exact shape of the filter function affected performance very little, 
the spread of the function had a large effect; the most suitable spread 
appears to be about three elements for the normal viewing distance. 

As the maximum permitted length of run is decreased from 10, it is 
found that there is very little increase in entropy for Code II for a 
maximum runlength even as short as 4, suggesting that a 4:1 grid 
algorithm may perform almost as well as free-running algorithms. 

The 2:1 grid algorithm (modified interpolative) does not permit 
operation at lower picture qualities and bit rates; the 4:1 algorithm 
has a larger range. However, over their range of operation, the grid 
algorithms perform at least as well as the best free-running algorithm 
and in view of their simpler implementation appear to be the most 
promising. 

Three different 2:1 grid algorithms were compared, an error-trans- 
mission technique in which the correction signal is sent as a difference 
between the estimate and the input, the modified interpolative al- 
gorithm, and the extrapolative algorithm. Extrapolation was slightly 
inferior to the other two methods and of these the modified interpola- 
tive method is more simply implemented. 

The emphasis in this paper has been on obtaining an efficient 
discrete representation of a picture signal rather than presenting a 
complete coding system. Consequently, there are a number of con- 
siderations such as sensitivity to transmission errors which are not 
discussed in the paper but nevertheless bear importantly on the feasi- 
bility of any practical coder. 
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Statistical Properties of Gilbert’s 
Burst Noise Model 


By MIN-TE CHAO 
(Manuscript received April 17, 1973) 


Simple statistical procedures for analyzing error data, e.g., in digital 
data transmission systems, are usually based on the assumption of in- 
dependence. This paper studies the performance and potential utility of 
such simple statistical procedures in the case of nonindependent error 
occurrences. The burst noise model is selected for this purpose because of 
ats neatness, its mathematical tractability, its built-in structure of de- 
pendence, and its importance in communication theory. We show that 
statistical procedures designed under the assumption of independence tend 
to be conservative for the burst noise model. For example, the usual bi- 
nomial test will reject, on the average, more channels with small error 
rates than it would if the errors were independent. The case that the sample 
size n and the error rate p converge in such a way that np — uo ts also 
studied. It is shown that the error process can be approximated by a 
compound Poisson process in continuous time t. The statistical implica- 
tions of this fact are also discussed. 


I, INTRODUCTION 


A dilemma long existing in the theory and applications of digital 
data transmission is the precise determination of the error structure. 
On the one hand, it is a well-recognized fact that errors do not occur 
independently; on the other hand, only the assumption of indepen- 
dence offers us a model sufficiently tractable that ordinary statistical 
procedures can be designed accordingly. A direct consequence is, 
of course, that we are using statistical methods designed for in- 
dependent observations to make statistical inferences on dependent 
data. 

The fact is, we do not have much knowledge of the error structure 
of data transmission channels. Mathematical models have been 
constructed for fitting observed data streams containing errors, 

18038 


1304. THE BELL SYSTEM TECHNICAL JOURNAL, OCTOBER 1973 


noticeably the burst noise model of Gilbert,! the Markov error process 
and renewal error process of Elliott,?? and the binary regenerative 
model of McCullough. 

One of the most pertinent models with a built-in dependence 
structure is Gilbert’s burst noise model. It is this model that we shall 
study in this paper. One of the prime concerns of this study is the 
behavior of various statistical procedures under the burst noise model. 

Gilbert! constructs a model for burst noise as follows. An input 
binary signal (0 or 1) is transmitted through a noisy channel with 
noise z (0 or 1) so that the output is given by 


output = input + 2 (mod 2). 


The channel can be in either of the two states, good (G) or bad (B). 
If, at time n, the channel is in G, there is no noise so z, = 0; if the 
channel is in B, a ‘‘coin’”’ with P[head] = h is tossed and z, = 1 is 
identified with a tail outcome. 

The channel can shift from a good state to a bad state and vice 
versa. Identify 1 as G and 2 as B and let X, denote that state of the 
channel at time n. It is assumed that the process {X,: 7 2 1} isa 
two-state Markov chain with stationary transition probabilities 


eR sl @ 


and initial distribution (7, 72). 

Let Z, = 21 +---+ 2, denote the number of errors through the 
nth-bit output (0 or 1) digits of the channel where z; = 1 if and only 
if an error occurs at the ith bit. The statistic Z, is obviously the 
quantity that will be used in any statistical procedure concerning the 
bit error rate. The statistical behavior of Z, will be studied extensively 
in this work. 

In Section II, we derive most of the exact formulas concerning Z,, 
including explicit expressions for its probability-generating function 
and its first and second moments. The exact form of the probability 
distribution of Z, is quite involved in general. For the special case 
p+P=1, Z, reduces to the binomial variable. The quantity A 
= 1 — p — P can thus be used as a measure of dependence; most of 
the complications in this work are caused by the presence of a nonzero 
\. The effect of dependence is discussed in some detail in Section III. 
Transmission in blocks of digits is considered ; one of our major results 
is that it can be shown in this model that the block error and the bit 
error have essentially the same covariance structure. Thus, most 
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results concerning bit error rate can be transferred easily to results 
about block error rate. As a corollary, the variance for Z, is obtained 
as a sum of two components, one due to the sum of variances (as if 
the z’s were independent) and the other due to the fact that \ + 0 
(the effect of dependence). 

Since Z, is known to be asymptotically normally distributed, the 
variance formula of Z, can be used to judge the effect of dependence 
on the robustness of statistical procedures (i.e., on how well procedures 
based on the independence assumption perform if this assumption is 
violated). A general conclusion of Section IV is that statistical pro- 
cedures designed under the assumption of independence tend to be 
conservative for the burst noise model. For example, the usual bi- 
nomial test will reject, on the average, more channels with small 
error rate then it is supposed to. 

It is shown in Section V that if the bit error rate p 0 in such a 
way that np — uo > 0, then Z, converges in distribution to a com- 
pound Poisson distribution. The statistical implications of this fact 
are also discussed. In particular, Z, is a minimal sufficient statistic for 
Lo(p) in some approximate sense. This justifies the use of Z, in any 
statistical decision procedures concerning the error rate p. 

Despite the model’s simplicity, the insight we gained in studying 
this burst noise model enables us to investigate more deeply the 
structure of error processes. For example, it is possible to treat the 
underlying Markov chain {X,} as an s-state stationary Markov chain. 
Details of this and other extensions and their implications will be 
discussed in a forthcoming report. 


II. STATISTICAL PROPERTIES OF Z, 


We shall assume, for simplicity and without loss of too much 
generality, in the sequel that the initial distribution (71, 72) of the 
two-state Markov chain {X,} agrees with its absolute stationary 
distribution (p/(p + P), P/(p + P)). Under this assumption, {X,} 
is strictly stationary. 

Let 

gn = P(Z, = 0]. 


Note that the bit error rate p is given by 


pS b= Hi 
= Pla = 1] 
(L = h)P/(p +P); (2) 
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and the block error rate p;, the probability that a block of size k 
contains at least one error, is 


pr = 1 — ge. (3) 
Thus, pi = p. 

Since the event [z; = 1] implies [X; = 2] and thus signifies a 
return to a bad state (a recurrent event), it is possible to utilize the 
renewal equation to derive an exact expression for g,. The following 
theorem is essentially due to Gilbert [Ref. 1, eq. (14) ]. 


Theorem 1: For n = 1, 
Ajapt! A sagt! 
1l—ai 1 — a?’ 








jr = (4) 


where 


a= elL— 1 = hy = Dp) = (ps = 2) 

+ VL — P) — AG — p)P + 4pPh] 
a2 = 3l— (1—A)(1 — p) — (p+ P — 2) 

=A (C= Py hao AapPr 
A, = play af (p eS 1) J/ar(ax = a2) 
Az = plaz + (p + P — 1) J/az(a2 — a1). 


A proof of Theorem 1 different from that of Gilbert (and the proofs 
of all other theorems) will be presented in the appendix. We remark 
here that since a broader view and a more systematic approach is 
adopted in our new proof, it is possible to extend our method readily to 
a more general framework than a two-state Markov chain. 

Relation (4) can be viewed as a relation between bit error rate and 
block error rate. If \ = 1 — p — P > 0, it can be shown that 0 < ap 
< a1 < 1 so that g, — 0 exponentially fast. One effect of dependence 
in this model is reflected in (4), namely that g, is the sum of two 
exponential terms instead of one. In general, if the underlying Markov 
chain is s-state, gn will be a sum of s exponential terms. 

The right-hand side of (4) is a function of p, P, and h. We shall 
write gn = gn(p, P, h) when we want to emphasize this point. An 
important connection between g, and Hu”, the probability-generating 
function (PGF) of Z,, is stated in Theorem 2. 


Theorem 2: The probability-generating function of Zn 1s given by 
Eu = g,(p, P, H), (5) 


where 


H= (1—A)ut+A. 
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Thus, replacing each h by (1 — h)u +h in (4), we obtain the PGF 
of Z,. The exact expressions for P[Z, = 7] are involved unless 7 is small. 
Using (4) and the fact that 0 < a2 < ai < 1, it is possible to express 
P[Z, = t] approximately in terms of its leading term as 





ee ee eo ee 
PLZ, = iw ft (") al . (6) 
Relation (6) can be used to establish the Poisson convergence of Zn 
if p = po/n > 0. However, an indirect proof will be presented later. 

Moments of Z, can be obtained by differentiating the right-hand 
side of (5) and setting wu = 1. Specifically, we have 


EZ, = Np, (7) 
7 = (n ~ Ud _ MC = A") 
Var Zn = no(1 — p) + 2c | a (a? I (8) 
where 
C = (1 = h)?ry12 
J=1-—-p-P. 


Relation (8) also can be obtained by oth»r methods which we shall 
discuss in Section III. 


III. MEASURE OF DEPENDENCE AND ITS EFFECT 


If the transition matrix of a Markov chain has identical rows, then 
this Markov chain is merely a sequence of independent and identically 
distributed (iid) random variables. For the two-state Markov chain 
{X,} underlying this burst noise model, the matrix 7 in (1) has 
identical rows if and only if p + P = 1. Letting \ = 1 — p — P, we 
see that |A| S$ 1 and that \ = 0 if and only if the channel is 
memoryless. 

The eigenvalues of the transition matrix play important roles in 
the theory of Markov chains. The largest (in absolute value) eigen- 
value is always 1; in general, it is the second largest eigenvalue that 
affects all the essential features of a Markov chain. The parameter \ 
defined earlier is the second largest eigenvalue of the matrix T in (1). 

The significance of the parameter \ can be interpreted intuitively. 
If p and P are small, the underlying Markov chain {.X,} tends to 
stay in a certain state (G or B) once it enters this state; hence, \ > 0 
indicates the tendency of producing burzsty errors. If both p and P 
are large, then {X,} tends to shift between the good state and bad 
state alternatively. Since the latter case is obviously not very interest- 
ing, we shall always assume \ = 0 in the sequel. 
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Let II denote the 2 X 2 matrix with identical rows 


71 12 
= ; 
TW1 72 


where (m1, 72) is the absolute stationary distribution of {X,}. By the 
definition of the absolute stationary distribution and by some simple 
calculations, it can be seen that 


UI? = TU = I? = I. (9) 
It follows from (9) and simple induction that, for n = 1, 
T* — Il = (T — I)” 


oes io 


Relation (10) allows us to calculate the f-step transition proba- 
bilities of {X,} accurately. It can also be used to find the covariance 
of z; and z;. We restate eq. (17) of Ref. 1 as follows: 


Theorem 8: The covariance of 2s, 2; (¢ *% j) is given by 


Cov (2;, 23) = CrAlH4, (11) 
where C _ (1 = h)?rime. 


Corollary: 
> = (os. le 
Define, for 7 = 1, 2, ---, 


T;=1 if 2 —1)k+1 + 2u-nyege Fees + Ze 2 1 
= 0 otherwise; (13) 


namely, 7; = 1 if and only if the 7th block of length & is not error-free. 
It is possible to extend eq. (11), and therefore (12), to the correspond- 
ing equations involving the T’s. 


Theorem 4: There exists 0 < Cy < © such that 
Cov (T;, Tj) = Cidl*-4*, (14) 


The value of C; can be found explicitly. However, we shall be satisfied 
with a crude estimate C1 = Corime\'—* where |C2| S t. 

Note that 7; = z; if k = 1. In this case, eq. (14) reduces to (11). 
Theorem 4 not only states that the 7T’s are ‘‘less dependent” than the 
z’s but it also tells us, in some sense, how much less dependent the 
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T’s are. Let 
Sn = Ti+ Totes Ds 


The statistic S, is the obvious statistical quantity to analyze if digits 
are transmitted in blocks of size k. For example, in the 1969-70 
Connection Survey*:® on the Bell System Switched Telecommuni- 
cations Network conducted by Bell Laboratories, statistics of block 
errors are presented for both high-speed and low-speed data trans- 
mission. Hence, the more important implication of Theorem 4 is that 
eq. (14) exhibits the same general structure as eq. (11). For example, 
replacing C by Ci, p by pz, and \ by A* = \* in (12), we immediately 
obtain the formula for Var (S,). 


Corollary: 
= * *2 __ \*n—1 
Var (Sn) = npx(1 Sale pr) als 204] (n 1)A = nN (1 ON ) 


hor To |’ me 


Consequently, statistical procedures using S, and concerning the 
inferences on the block error rate p; should have essentially the same 
behavior as those procedures using Z, and concerning the bit error 
rate p. The above reasoning implies that, at least as far as the large 
sample properties are concerned, it is sufficient to consider inference 
on p only. 

Both the law of large numbers and the central limit theorem hold 
true for the sum of Markovian random variables; see, for example, 
Ref. 7. Hence, 


Sa 
Gereee (16) 
with probability 1; and 
Sn — Npz | 
P| ———— £v|—- ®(v 17 
WEES: ay ue 
for each — «7 <v < o, where @(v) denotes the cumulative distri- 


bution function of an N(0,1) random variable. Relations (16) and 
(17) will be used in Section IV to discuss the robustness of some 
statistical procedures concerning inferences on px. 


IV. STATISTICAL INFERENCES ON px 


For simplicity, we shall consider the special case k = 1 and concen- 
trate our discussion on problems of statistical estimation and hy- 
pothesis testing of p = pi. As remarked earlier in Section III, the 
restriction k = 1 can easily be extended to the general case. 
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Since {X,} is assumed to be stationary, so is {zn}; we have seen that 
ElZ,] = np, (18) 


so that the obvious estimator 6, = Z,/n of p is unbiased. Relation 
(16), specialized for the case k = 1, states that 6, is a strongly con- 
sistent estimator of p. 

Very few (optimal) small sample properties of 6, can be stated, 
however. For n 2 38, it can be shown that no uniformly minimum 
variance estimator of p exists. Nevertheless, it is intuitively obvious 
that fn is about the best we can do if the 2’s are the only observables. 
From (12), 


n Var (Bx) = p(1 — p) + 207~~ + 0(1) (19) 
4 e+ A, 
where 
o? = p(1 — p) 
A= 2(1 _ h)*ayred/ (1 se d). 


Note that the term p(1 — p) in eq. (19) corresponds to n Var (f,) if 
the z’s were independent. Since we have assumed that \ = 0, it follows 
that A = 0 and n Var (62) 2 o?. Thus, the presence of a positive 
actually causes loss of efficiency in estimating p. Writing A = 7’, 
we see that if the parameters h, p, and P (hence o? and 7?) can be 
estimated from the data, the loss of efficiency due to dependence can 
be estimated as the ratio #/¢, where 7 and ¢ denote the estimates of 
7 and o from the sample. Hence, if control or confidence limits are 
used to evaluate the channel performance, the actual 3 standard 
deviation or 2 standard deviation limits should be wider by 100(7/c) 
percent. 

We may also consider the loss of power for statistical tests for Ho: 
p & po of the form 

reject Hp if Z, 2 C*. 


Based on the assumption of independence, the power function is 


approximately 
c* — ve ) 
=] —o( —__* }, 20 
7 (es oO 


whereas for our model, the power is approximately 


6 = 1— 2( Tee). (21) 
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P[Z,2C*|p] 





BIT ERROR RATE p 


Fig. 1—Comparison of power functions. 


If the first-type error a S$ 34, we see from (20) that C* — npo 2 0. We 


see that Br S Bp if C* — np 2 O and £B; 2 Bp otherwise. This means 
that it might be possible to design more powerful tests for Hy based 
on the knowledge that the dependent model obtains. On the other 
hand, the test is conservative in the sense that it may reject more 
channels than expected if the bit error rate p is close to the service 
objective po and if the dependent model obtains. The rules of the game 
shift in the other direction if C* — np < 0. However, it is the smaller 
values of p that we are really concerned with and we may claim that 
the test based on the assumption of independence gives a pessimistic 
estimate of channel reliability (see Fig. 1). 


V. POISSON APPROXIMATIONS 


The bit error rates of high-speed digital channels are usually very 
small, say 10-5 to 10-8; therefore, the normal approximation and the 
statistical theory discussed earlier may not be too helpful in practice 
unless n is large. In this section, we prove that Z, converges in distri- 
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bution to a Poisson distribution if np — yo in a suitable way. Using 
this result, we construct a Poisson process in continuous time ¢ that 
approximates the process {Z, (¢):t > 0} where n denotes the number of 
transmitted digits per unit time. 

We have shown earlier in (2) that the error rate p is given by 


pal Pipe P): (22) 


If p —0 in such a way that np — uo > 0, what do we expect to be the 
asymptotic distribution for Z, = 21 +---+ 2n, the number of errors 
in the first n digits? Note that we have quite a few choices for the 
convergence np — uo. For example, keeping p fixed and letting 
P= (u/n)",1 —h = (u/n)®, €: + eg = 1, we have, by (8), 





Var (Zn) & wo(1 — 9) + 2Cn-s us . 
2 1- 
uot Satta, (23) 
where uo = u/p. Also, 
EZn = np 
— Lo. 


Hence, if eg = 0 is selected, we see that for large n, Var Zn ¥ EZn 
so that the limiting distribution of Z, cannot be Poisson. 

In order that p = (1 — h)P/(p +P) Suo/n, the most general 
choice of h and P would be 


1 —h = aye + aox? + aga? +---, 
P = by + boy? + bsy? +:--, 


where « = n-®%, y = n-%, ey + eg = 1, 1: 2 0, €2 > O, and aybi/p = po 


(the case e2 = 0 is of particular interest and will be considered sepa- 
rately later). We state the main theorem of this section as follows: 


(24) 


Theorem 5: If p—0 in such a way that (24) holds, then 
P({Z,=t|]- 5 poe *° 


as n— ©, where uo = Qibi/p. Furthermore, the convergence is uniform 
int = 0,1, 2, -->. 

By using the result of Theorem 5, we may construct a Poisson 
process in continuous time ¢ as an approximation to the process of 
partial sums {Z,: » 2 1}. Suppose the underlying channel can 
transmit n digits per unit time. Let Z,(¢) denote the number of errors 
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in (0, t). Theorem 5 states that, for 7 = 0, 1, 2, ---, 
PLZ a(t) = tJ F (ule; 


here » denotes the limiting error rate per unit time. Let Z(t) denote 
the number of errors in (0, t) in the limiting case. The fact that Z(t) 
is a process with independent increments, namely that Z(t) is indeed 
a Poisson process, is easy to prove and we shall omit it. 

Theorem 5 implies that Z, is asymptotically a minimum sufficient 
statistic for the bit error rate p if (24) can be justified; this provides 
theoretical support for the use of Z, in any statistical inferences 
concerning p. We remark here that, by replacing p by p; and Z, by Sp, 
the same comment applies for block error rates. Another consequence 
of Theorem 5 is that 

Var (Zn) — Lo = aybi/p. (25) 


Note that if \ = 1 — p — P = O (the independent case), (24) implies 
P — 0 and this in turn implies p — 1. From (25), we see that Var (Z,) 
is minimized in the independent case. The increase of variance due to 
dependence is therefore 100(1/p — 1) percent. Hence, in the dependent 
case, the confidence interval for p should be wider than we thought 
in the independent case. 

The null hypothesis Ho: p S po becomes Ho: uo S uo in the limiting 
case. The uniformly most powerful test for Ho exists and is given by 
the rule: 


reject Hy if Z, = C*. 


Based on the approximation that Z, is Poisson, we may compute the 
power functions as 


Bp (uo) 


I 


P[Z, = C*| uo] 


re 
= — ure #0 
pes a! Bo 


zis f — 2, c*—1q 
= , (Tp! D1 ee x 
and 


= as 1 — ty C*—1d 
p= f Chat” © B. 


It follows that 8; S Bp so that a test for Hy based on the assumption 
of independence and used when dependence is present rejects more 
channels than it should. In other words, tests designed for independent 
observations protect customers in the sense that channels they are 


using may have better quality than inferred. 
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The effect of dependence reported for both the binomial and the 
Poisson cases has an intuitive explanation. By using Z, or S,, we are 
actually abandoning some of the information contained in the sequence 
Z1, 22, -*-, So that statistical inferences based on Z, or S, tend to be 
more conservative in the sense that channel reliability is estimated 
pessimistically. 

We now return to (24) and consider the special case ez = 0. This 
case cannot be ignored because previous papers, for example Ref. 1, 
indicate that sometimes h 0.5 (rather than 0.999) is a reasonable 
value. The fact that Poisson processes do not describe certain error 
processes well has also been reported in the literature. 

If eg = 0, eq. (24) reduces to 


P = [61 + 0(1) J/n bi > 0. (26) 
We have 
Theorem 6: If (26) holds, then 
bil — H) 
iq, Zn —{—] 
Eu — exp | (rae | (27) 


where H = (1—h)u +h. 

We remark here that the limiting value in eq. (27) is the PGF of a 
compound Poisson process. More specifically, let N be a Poisson 
variable with mean 61/(1 — p), and let Wi, We, --- be iid random 
variables with the geometric distribution 


7=0,1,2,---. 

If the W’s are independent of N, then the left-hand side of (27) is 
simply Hu¥1+W2+---+Wv, It is of course possible to introduce a con- 
tinuous time parameter ¢ and consider the following random mecha- 
nism which describes the bursty nature of this error process vividly. 
The bursts are generated by a Poisson process; given that a burst 
occurs, the errors are generated by a geometric distribution. 

From the right-hand side of (27), it is possible to compute the 
moments of the limiting distribution of Z,. We have 


bi(1 — h) 


- (29) 


E(Wi+ W2+--:+ Wy) = 
Var (Wi + We.+::-+ Wy) 
a bi(1 a h) te 2bi(1 = Ot = Dp). (30) 
Pp P 
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Note that the variance is always larger than the mean in this case. 
Note also that, as h approaches 1, the second term on the right-hand 
side of eq. (80) is of higher order and vanishes in the limiting case. 
Another interesting thing is that it is possible to show that the right- 
hand side of (30) is minimized at p = 1, and as p approaches 1, the 
limiting distribution is Poisson. 

Branching renewal processes have been suggested in the literature? 
as a model for series of events. The basic structure for branching 
renewal processes can be described in terms of our problem as follows: 
The series of primary events (bursts) are generated by a Poisson 
process. Each of these primary events generates a subsidiary series 
of events (bit errors), separated by the waiting time Y,, Y2, ---, Ys, 
where S is random. If we assume that these subsidiary series of events 
take no time, then the branching renewal process reduces to the com- 
pound Poisson process. 


VI. CONCLUSIONS AND FURTHER EXTENSIONS 


(1) The burst noise model of Gilbert discussed in this paper pro- 
vides a vehicle for studying the robustness of some fixed sample size 
statistical procedures. The general result is that the presence of 
dependence increases the variance of the random variable Z,, for 
the case where the bit error rate p is fixed and the case in which p 
= [uo + o(1) ]/n. Thus, use of statistical tests based on the assump- 
tion of independence increases the power at the cost of rejecting more 
satisfactory channels than would be rejected if dependence were 
absent. The use of blocks does reduce the covariances among errors 
compared with bits or smaller blocks. However, the covariance 
structure among the blocks is essentially the same as that among 
the bits. 

(77) Although the dependence structure of the Gilbert’s burst noise 
model is a simple one, it is by no means a trivial one. In fact, from the 
insight gained through this study, many results obtained in this paper 
have generalizations in error processes defined over an s-state Markov 
chain as well. A unified treatment on channels with Markov type of 
memory will be reported elsewhere. 

(tat) The second largest eigenvalue (in absolute value) of the 
(s X s) transition matrix of the underlying Markov chain is a param- 
eter which should not be overlooked. It can be viewed as a measure 
of dependence of a Markovian model. The effect of this parameter 
(=) in this work) is visible in many important formulas, for example, 
in (14). 
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(cv) Another important question to ask is what kind of stochastic 
process can be used to approximate the error process of a binary 
channel with memory. If the bit error rate is small, we can extend the 
proof of Theorem 6 (in a nontrivial way) to find an important con- 
clusion: the compound Poisson process can serve the purpose. 

(v) The by-products of this work are also fruitful. For example, 
the variance formula of Z, can be generalized to find the variance of 
Tn = f(X1) +-+--+ f(X,) where {X;} is an s-state Markov chain, 
s S o, and fis an arbitrary function. Since many continuous sampling 
plans, such as CSP1, CSP2, CSP3, can be described as random walks 
of the form 7, (see Refs. 9 and 10), the application of this formula to 
quality assurance is evident. 

(vt) Mathematically speaking, there is an essential difference be- 
tween Gilbert’s original treatment and our generalizations to the 
s-state Markov chain. More specifically, Gilbert viewed his problem 
as one of the renewal type whereas the s-state Markov case should 
be handled by the semigroup property (of taboo probabilities). We 
remark here that many results of the theory of recurrent events (see, 
for example, Ref. 11) can be applied to Gilbert’s model. We also 
remark that the renewal process is a one-state semi-Markov process. 
A general question can be raised at this point: What is the behavior 
of an s-state semi-Markov channel? Since it is known that distri- 
butions other than the exponential (for example, the Pareto distri- 
bution, see Ref. 12) describe the waiting time distribution well, the 
question raised is a realistic one and should not be merely considered 
as an attempt at mathematical generality. 
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APPENDIX 
A.1 Proof of Theorem 1 


Consider Y, = (Xn, Zn) as a three-state Markov chain with transi- 
tion matrix 


(G,0) —_(B, 0) (B, 1) 
(G, 0) 1-—P hP (1 — h)P 
(B, 0) p hAl—p) (l-AA—p)|=2= (as), (81) 
(B, 1) p hep) “Ch =f) =p) 


say. We have 
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Buy4n = Hyttteets ten 


= sare ee is sees De gaan EEO ysis 
aa! Paar a ae Qyour' yo» (32) 


where \,, = PL Yo = yo]. Note that the value of z; is completely de- 
termined by the value Y, = y,. Let 


Py pay = Udy uy 


and let 
R = (753). 
Relation (32) can then be written as 
Eu*® = VR"1, (33) 
where 
WM = (AG,0), 4B,0), AB), 
P= ,1.4), 
R = (rij) 
1 0 0 
=@Q@10 1 O}- 
0 0 wu 


We remark here that eq. (33) can be extended to the case of an s-state 
Markov chain easily. Letting u = 0 in eq. (33), the PGF of Z,, we have 


1 0 0\} 
P[Z, =0] = (s 1 7 i (34) 
00 0 


The last column of the 3 X 3 matrix in eq. (34) is always a zero vector 
for every n = 1. Hence, the right-hand side of (34) is essentially the 
nth power of a 2 X 2 matrix. The explicit formula for g, in eq. (4) 
follows from (34) by straightforward calculations. 


A.2 Proof of Theorem 2 


The z’s are conditionally independent if the values of the X’s are 
given. Hence, 


PLZ, = 0] = E{ Plz: = 2 = +++ = 2, = 0|X1, Xe, vos, Xa} 
= BUTT Pla = 0/3] 
es 


= BY Wk 
t=1 
= A{pXi1t Xo 4X pon, (35) 
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Similarly, 
Ey = E{ E[uztet- eben |X4, Xo, ree, Xn]} 


Il BLw XJ} 


P| 
BlT eae (i= muy 
i= 
= BYXA Kat Xamon (36) 


where H = h + (1 — h)u. By comparing eqs. (85) and (36), Theorem 
2 follows. 


A.3 Proof of Theorem 3 
By eq. (10), 
P(X, = 2[|Xo = 2] = mo + Am. 
Hence, 
Cov (z:, 2;) = Plz: = 2; = 1] —- ?’? 

= (1 — A)?PLX; = X; = 2] — 
(1 — h)?2o(1e + wyA!*-J1) — p? 
mimo(1 — h)?A!*41, QED. 


A.4 Proof of Theorem 4 


Let us compute a special case first. Consider P[T, = 0, 7, = 0]. 
A typical path of the underlying Markov chain {X1, Xo, ---, Xin} 
may be of the following form: 


D(X, Capt, Ln) ky L(n-1) k+1) 
= (1%2°° URC REL? *§ ° XL (n—1) kL (n—1) k4- 10 (n—1) 42° ° *“2nk)- (37) 
| 


first block last block 
of size k of size k 


The rest of the x’s in b (and in Wi, W2 later) are omitted for typo- 
graphical reasons. Note that the values of 2x42, ++, Z(n—1)r-1 are 
deliberately unspecified ; also, n > 2 is assumed pro tem. 

For fixed first block and last block, there are four different kinds of 
paths, according to the values of #441, %(n—1)x. Let m denote the number 
of 2’s in the first and last blocks together. We have 


P[T, = 0, Ta = O[D(tz, Levi, Vom—1) ky Lin—1) e411) |] = A™ (38) 
and 


PLB (ee, Seer, Lon=1)% Lea) | 
= Wi (CE) PaesPiee i Vee apeeaen 4 (2¢n—1) k+ 1); (39) 
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where 


W(x) = PLX, = Vy, Xo = La, **"y Xjp-1 = Vk-1; Xt = xz | 
and 
Wo(a(n—1) +1) = PLEX (n—1) k42 = L(n—-1)k+2) °° %y Ng = a 
|X (n—1) k41 = Bin rye. 

We may also find P[7T, = 0]-P[T, = 0] by considering their 
conditional probabilities over the first and the nth blocks. It is not 
difficult to see that 

PLT; = OJPLT. = 0] = 23 h”W 1 (xx) Wo(%(n—1) b+) * Pa¢n_ayegas (40) 


where the summation ranges over all 2?* possible blocks. The expression 
for P[T: = 0,7, = 0] can be obtained by taking the product of 
(88) and (39) and summing over all 2?*+? possibilities. The 274+? terms 
in this form of P[7T; = 0, T, = 0] outnumbered the terms in (40) 
by a margin of 4 to 1, and there is an obvious 4:1 correspondence 
between these terms. Consider 
Cov (T1, Tn) = Cov (1 = T 1, {= Tn) 
= P[T,=0,7, =0]—P[T, =0]P([T,=0]. (41) 
For fixed first and last blocks, a typical difference between the (4:1) 
correspondent terms is 
hm Wi1(xx) Walter) pep Opie iia 
PRs ap PD aeiactnas + Depa Deen 
— Dospan Dasieciga— ie eesyeenle (42) 


By (10), it can be shown that the third factor of (42) becomes 


—2) k—2 
(D2,1P12¢n-1)k4t a Px, lP22n—1yk41 a Pxy2P 12 m_ry E41 + D2, 2P22¢4-1ynu1) 0 
P 





poe epee Sak (ty Zana) = 1, D) 

= et  (n—2) k+1 = (1, 2) 

= — ) (n—2) kt = (2,1) oe) 
Pp 
; 7 3 \(n-2) +1 = (2, 2). 


Note that in all terms we have a common factor \‘"-®)*+!. By factoring 
out this common factor, Theorem 4 follows immediately. 
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We may even push the computations further to find an exact 
expression for the constant C; in Theorem 4. Note that we have four 
types of combinations of blocks, according to the values of x, and 
XL(n—k)+1. The quantity in (42) becomes 


P 
(1, 1) => h"Wi(1)W2(1) oTP \ (m2) HL, 











(1, 2) => — he W3(1) We(2) —— yore, 
pede 
(44) 
—. fm (n—2) k+1 
(2, Le hmW1(2)W2(1) D+P nN ; 
2, 2) => h™W(2)W2(2) Pret 
(2, 2) => mW (2) Wa(2) Pg none, 
and Cov (71, T',) is the sum of all 2?* terms in (44). 


Let 
Ty = #1 + 22 bee + en, T. = 2np2 tees Zon. 


(We should use T, = 2(n—-1yez2 +:°+++ Znzj; however, the distribution 
of T,, is independent of n so we may take n = 2.) Then 


PCT; = O|X;, = LIPETS = O|Xi+1 = 1] 
a hm P(X, = 21, vee Xpy = Xz1| Xp = 1] 


22k-1) 


*PUX nyo = Lape, +++, Xan = Lox|Xeyi = 1] 


1 
= DS he W yA) Wad); 
WY Q2k-1) 
where m denotes the number of 2’s in the sequence 21, %2, -*+, 24-1, 
Lk+2, ***, Cox, Which is equal to the number of 2’s in the sequence 1, 22, 


+++, 22, in the case %, = 241 = 1. Thus, the sum of terms of the type 
(1, 1) in (44) is simply 


mP[T, = 0|Xz = LIPLT2 = 0[ Xia = 1] -merorver, 


Similarly, we may find the sums of other types of terms in (44). 
We have 


(1,2) => — mAP[T) = O|X, = 1]P[T2 = 0|Xuga = Q]rod(r-2) eH 
(2,1) = — rehP[T, = 0|X, = 2]P[T2 = O[ Xanga = 1] 41 
(2,2) => weh?P[T, = O|X_, = 2IP[T2 = O| Xn = Q]rir-*4, 


I 
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Thus, ifn > 2,7 = 1, 


Cov (T1, Tn) 
= Cov (T:, T it-n—1) 
= mired" PTT, = O|X, = 1LIP[T: = 0/ Xai = 1] 
- hP[T, = 0|X, = 1)P(T2 = 0/ Xin = 2] 
= hPLT; = 0|X; = 2)P(T» = 0/Xi1 = 1] 
+ WPCT, = 0|X, = 2)P[T. = 0! Xas1 = 2]}. (45) 
The case n = 2 should be considered separately; this is because 
(n — 2)k —1< Oif n = 280 that (39) simplifies to 


PLb(xx, 41) ] = Wil&e) Pepeps: WV 2(ee41). (46) 


In this case, the number of terms in P[T: = 0, Tz = 0] equals the 
number of terms in the product P[7T:1 = 0]P[T: = 0] and there is 
an obvious one-to-one correspondence between the terms. Consider the 
difference P[T, = 0, T2 = 0] — P[T: = 0] PLT. = 0]. For fixed 
first and last (second) blocks, the term-wise difference is 











h™W(xx)We(enr1)[Dopeng: — Tred (47) 
The last factor in (47) can be computed. We have 
Popena — Tey = iP hf (ay Tapa) = (1, 1) 
=-5 5) = (1,2) 
aia + 5 = (2,1) 
= 5 F = (2,2). (48) 


Note that (43) reduces to (48) if n = 2; hence, all arguments leading 
to (45) hold true even if n = 2. 

Let C2 be the quantity in the large square bracket of (45); we may 

write (45) as 
Cov (T:, T;) => Comyrod (lt s ITD) 1 (49) 
aS. 

It is possible to find the value of C2 through an argument similar to 
that of finding g, in Theorem 1. However, we shall be satisfied with a 
crude estimate 

| wyreC's | < Lis 


which follows from mir2 = 71(1 — 71) S 2? trivially. 
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A.5 Proof of Theorem 6 


Let H = (1 — Aju +h and let 41, &, A1, As, p be the quantities 
obtained from a1, a2, Ai, A2, p by replacing each h with H respectively. 
Asn—o, H —1. Hence, 

ai 1 
Ax — 0 
As — 0 
p— 0. 
It follows from Theorem 2 that 
fA n+1 
inne toe (50) 
n> 0 no — GQ 

Let A? = [1 — P — H(1 — p) }?? + 4pPA in the expression of @1. 
An important step in our argument is to find the value of A. By 
substituting (24) into the expression for A*, we have 





a= pit Y yet + Y dy" + ey + olay), (61) 


n= 


Yon = (1 — p)?(1 — u)?(an + 2aiden—1 + 2d2G2n-2 ++ °° 
+ 2dn1An41) + 2p(1 — p)(1 — u)aen 
Yang = (1 — p)?(1 — uw)? (Qaiden + 2a2den1 ++ 2Andn41) 
Son = D5 + Qiben—1 + Zbaban—2 +++ 2bn-10n41 
bongt = 2dydon + Zobon-1 Feet 2bnbngi 
e= — 2(1 — p)(1 — wards — 4p(1 — w)arbr. 


A=pt X (dau® + eny") + fay + o(zy). (52) 
By comparing the A? in (52) with the same quantity in (51), it is not 
difficult to see that 


d, = (1 — p)(1 — u)ax 
ey = bi 


(53) 
for k = 1, 2, ---. Also, it is easy to find that 
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Using (53), it can be seen that in the expression of &:, the coefficients 
of x*, y* are zero for all k. Hence, recall that xy = 1/n, 


&=1- ewan xy + o(xy) 
ee ee come LE (=): (55) 
np n 
By (55) and (24), it can be seen that 
4,- 02d) 1 o( =): (56) 
pn n 


By (50), (55), and (56), we have 


lim Bu?" = exp (ee= 2) 


which is the PGF of the Poisson distribution with mean equal to 
ayb,/ D. : 


A.6 Proof of Theorem 6 


Using (26), we may express @1, @1 in terms of powers of 1/n as 





jy ee I 1 
ea = +0(:), (57) 
Ay = bes + o( =) ; 
n n 
where 
= bi =A) 
“" T—H-+ oH 98) 
H=(1-—A)uth. 
By eqs. (50), (57), and (58), we have 
litt But = lim —“2= lim ay" 
no n> 0 1] — Ai noo 
= exp[—a] 
- oe DUC Seen) oe 
= exp] — re pe ae 


This proves Theorem 6. 
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The Metallurgy of Remendur: 


Effects of Processing Variations 


By M. R. PINNEL and J. E. BENNETT 
(Manuscript received May 14, 1973) 


A recent development effort in telecommunications switching apparatus 
has been directed toward the production of a remanent reed, dry, sealed 
contact (remreed). Remendur, a medium-hard magnetic alloy nominally 
composed of equal parts tron and cobalt and 2.7-wt. percent vanadium, 
was chosen as the reed material in this contact. However, the application 
required the alloy to possess rather specific magnetic and mechanical 
properties and considerable difficulty was experienced in consistently 
processing Remendur into wire with these specified properties. To as- 
certain the sensitivity of these properties to variations in processing times 
and temperatures, and vanadium content, two melts of Remendur (2.5- 
percent V and 3.0-percent V) were processed with selected alterations in 
annealing temperatures at several stages. Microstructures were char- 
acterized following each step by light microscopy and were correlated with 
the appropriate ternary equilibrium diagram. Results demonstrate that 
microstructures developed by anneals between 900°C and 950°C are 
extremely sensitive to the precise temperature of the anneal and compost- 
tion of the alloy. The microstructure, which strongly influences magnetic 
and mechanical properties, can be varied over the limits of the two-phase 
ai + y region by variations in vanadium content of only 0.5 wt. percent 
and by the small 50°C temperature range. 


I. INTRODUCTION 


Historically, considerable difficulty has been experienced in con- 
sistently processing the Fe-Co-2.7-wt. percent V alloy (Remendur) 
into wire with specified magnetic and mechanical properties. A major 
cause for these problems has been that the metallurgy of this ternary 
system was not sufficiently understood. A number of investigations 
have provided information on various aspects of phase equilibria!~* and 
kinetics and mechanisms of phase transformations®-” in this system. 
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The most notable of these are the work of Késter and Schmid! on the 
phase equilibrium of the system and the work of English® and Martin 
and Geisler? on the (FCC) S (BCC) transformation. However, dis- 
crepancies still exist on aspects of phase equilibrium,!:?+;7 mechanical 
ductility,?-°7-" transformation kinetics,?:*7° and the influence of these 
factors on the development of magnetic properties.?:*7:19.!2.13 These 
discrepancies made analysis of the sensitivity of the mechanical and 
magnetic properties to composition and to heat treating times and 
temperatures during processing unreliable with further investigation. 
This analysis is vital since 2.7-percent V Remendur is being used as 
the reed material in the remreed contact. This paper describes the 
characterization of low-vanadium (2-3-wt. percent) Remendur, mainly 
by light microscopy, for the various processing steps from 6.35-mm 
(0.25-inch) rod through 0.538-mm (0.021-inch) wire to 0.18-mm 
(0.007-inch) flattened strip. Correlation of the microstructures with 
the equilibrium phase diagram of Késter and Schmid! is provided. 
Aspects of the microstructures in relation to cold ductility and mag- 
netic properties are discussed. 


II. MATERIALS AND EXPERIMENTAL PROCEDURES 


Two melts of Remendur, each with a different vanadium content, 
were used in these experiments. The alloys were prepared by Battelle 
Memorial Institute with nominal compositions of 2.5-percent V-balance 
equal Fe/Co and 3-percent V-balance equal Fe/Co. The actual 
analyses as determined by atomic absorption spectroscopy are given 
in Table I. The typical Remendur processing sequence from melting 
through final reed fabrication is outlined in Fig. 1. The influence on 
cold ductility, yield and tensile strength, resistivity, and magnetic 
properties of several steps in this processing sequence was unclear. 
These steps included temperature of and rate of cooling from the 
6.35-mm rod anneal, the need for intermediate strand anneals, and the 
temperature of the 0.53-mm wire strand anneal. Also, the influence of 
the stamping operation relative to the unworked shank on final reed 
properties was of interest. 


TABLE I—REMENDUR COMPOSITIONAL ANALYSES 
(wt. %) 
V Co Fe 


3% V Battelle 2.97 48.70 47.34 
2.5% V Battelle 2.46 48.36 48.27 
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6.35 mm ROD ANNEAL 


MELTING — CASTING HOT ROLL TO 6.35 mm ROD 
Fe — 48.5% Co —2.7% V 1100 — 1200°C IN AIR, Bae etnies: eRe ES 
3x3 INCH INGOTS AIR COOLED IN ICED BRINE 
SURFACE CONDITIONING 
OF 6.35 mm ROD COLD DRAWN TO INTERMEDIATE STRAND 
MECHANICAL AND/OR 1.27 — 1.52 mm WIRE ANNEALS (OPTIONAL) 
CHEMICAL 
COLD FINAL STRAND ANNEAL STAMP AGING ANNEAL 
DRAWN 900 — 950°C IN DRY Hp, AND 600 — 620°C IN DRY 
TO 0.53 mm 30 SECONDS, COOLED IN CUT Ho,2 HOURS, FURNACE 
WIR WATER-CHILLED H, CHAMBER REEDS COOLED IN H, 


Fig. 1—Typical Remendur processing. 


An experimental program was developed to provide clarification of 
the influences of these factors. A synopsis of this program is given in 
Fig. 2. It involved carrying two compositions of Remendur through the 
entire processing sequence with appropriate variations in parameters 
at each critical stage. For temperature variation, a matrix of tem- 
peratures between 900°C and 950°C was investigated. This range was 
based on previous statements by Gould and Wenny that cold ductility 
could be obtained by an iced brine quench from 925°C” and the current 
use of this temperature by K. Olsen (Bell Laboratories-Murray Hill) 
for wire processing.1® 

Microstructures were evaluated primarily by optical microscopy. 
Metallographic preparation was routine with a 5-percent Nital solu- 
tion used to lightly etch (10-50 seconds) the polished surfaces. Ob- 
servation and photography were carried out on a Ziess Ultraphot III 
metallograph using Nomarski Differential Interference Contrast 
(DIC). Indications of the ductility of both hot-rolled and annealed 
6.35-mm rod were provided by the drawability of the material. 


III]. RESULTS AND DISCUSSION 
3.1. Hot-Rolled 6.85-mm (0.25-inch) Rod 


A vertical section through a portion of the Fe-Co-V ternary equilib- 
rium phase diagram (Késter and Schmid)! near the nominal com- 
position of Remendur is given in Fig. 3. At 2.7-percent vanadium and 
temperatures above 950°C the alloy is single-phase FCC (vy). This is 
the structure of Remendur during hot rolling. The typical micro- 
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2.5% V ALLOY 3.0% V ALLOY 


1200°C HOT ROLLED 6.35 mm ROD 1200 °C 


— MICROSTRUCTURE 
— INFLUENCE OF COOLING MODE 
— DRAWABILITY 


900° C 925°C 950°C ANNEALED 6.35 mm ROD 900 °C 925°C 950°C 


— MICROSTRUCTURE 
— DRAWABILITY 


900 950900 950 900 950 900 950900 950900 950 


925 925 925 COLD DRAWN AND STRAND 925 925 925 
ANNEALED 0.53 mm WIRE 


— MICROSTRUCTURE 

— RESISTIVITY 

— MAGNETIC PROPERTIES 

— MECHANICAL PROPERTIES 


600°C AGING ANNEAL OF ROUND 600 °C 
AND FLATTENED 0.53 mm WIRE 


MICROSTRUCTURE 
RESISTIVITY 

MAGNETIC PROPERTIES 
MECHANICAL PROPERTIES 


Fig. 2—Synopsis of Remendur processing experiment. 


structure developed by cooling from the hot-rolling stage is shown in 
Fig. 5a. This microstructure is irregular and unconventional. English® 
and Chen’7'8 have amply reviewed the characteristics of this structure 
and only the most pertinent factors will be restated. Although this com- 
position passes into the two-phase a; + vy field on cooling, the trans- 
formation is sluggish. As a result, with rapid cooling a nonequilibrium 
single-phase structure which is entirely BCC (a1), supersaturated in 
vanadium, is developed (Fig. 5a). No retained y was detected by 
X-ray analysis. This has been referred to as a “‘massive”’ transforma- 
tion by English® and designated as the a2 structure. This phase, being 
a nonequilibrium structure, does not appear on a conventional equilib- 
rium phase diagram. Chen7® proposed that the transformation is 
actually the more common “martensitic” type. However, the structure 
of the transformed material is similar for both types® and, hence, this 
will not be discussed here. 

The ductility of Remendur with the all a2 microstructure and as air 
cooled was tested in a qualitative manner. An attempt was made to 
draw 6.35-mm rod in this condition (i.e., as hot rolled) into 0.53-mm 
wire by successive reductions of approximately 20 percent in area. 
The rod could not survive the first draw pass due to extreme brittleness. 
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Fig. 3—Vertical section through Fe-Co-V ternary phase diagram, from the work 
of Késter and Schmid. 


3.2 Annealed 6.385-mm (0.25-inch) Rod 


The brittleness of the air-cooled a2 microstructure dictates the need 
for an additional heat treatment to create drawability. A treatment of 
one-half hour at 925°C (in the two-phase a; + vy field) followed by an 
iced brine quench had been found satisfactory in producing cold 
ductility.!2-15 However, the sensitivity of ductility to this precise time 
and temperature and the influence of the structure produced by this 
anneal on subsequent properties developed in the drawn wire had 
yet to be investigated. 

Samples of both the 2.5-percent V and 3.0-percent V alloys in the 
hot-rolled and air-cooled condition were annealed at 950°C, 925°C, 
and 900°C for one-half hour and quenched in iced brine. The cor- 
responding microstructures are shown in Figs. 4 and 5. It is apparent 
that a significant variation in structure occurs as a function of tem- 
perature over the relatively small range of 900 to 950°C. At the same 
annealing temperature, a significant variation as a function of the 
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Fig. 4—Quenched microstructures of 2.5% V-6.35-mm Remendur rod in hot- 
rolled and various annealed conditions; etched; DIC; 600X. (a) Hot-rolled. (b) 
eee ; $ hour; 950°C. (c) Annealed; 4 hour; 925°C. (d) Annealed; 3 hour; 


small difference in vanadium composition (2.5 to 3.0 percent) is also 
obvious. 

Analysis of the equilibrium phase diagram for this ternary system is 
necessary to understand the critical nature of temperature and com- 
position for transformations in this range. Based on the vertical sec- 
tions of the ternary space diagram provided by Késter and Schmid,} 
isotherms (horizontal sections) at temperatures of interest were con- 
structed by the authors. The Fe-Co-rich portions of isotherms at 950, 
925, and 900°C are drawn in Fig. 6. The nominal composition of 
Remendur (2.7-wt. percent V-balance equal Fe/Co) is represented by 
an X. At all three temperatures the equilibrium structure is in the 
two-phase (a; + y) region, but the amount of each phase can be seen 
to change radically with temperature. This is a consequence of the 
y phase field shrinking rapidly toward the Co corner and the ai + y 
phase field increasing significantly in area as the temperature decreases 


PROCESSING OF REMENDUR 1331 








Be oe at ‘ Aen aa, a 


Pinay 
ey 3) Sh ey 

fh \% » 
ey a ~e ‘ 
q 4 pon SEEN 


aug 
<\aSSeES ae 
ny Stax omameaes 











wey SN SS = SAR: ne Sh re 
YAI 2 3 foe 

\ at Pies ROS SS ‘ Ne 
gi Ys oy ) SaaS eta Be 4 At 


et head d Po 





Fig. 5—Quenched microstructures of 3.0% V-6.35-mm Remendur rod in hot- 
rolled and various annealed conditions; etched; DIC; 600X. (a) Hot-rolled. (b) 
coe 4 hour; 950°C. (c) Annealed; 4 hour; 925°C. (d) Annealed; } hour; 


from 950°C to 900°C. This a; phase (BCC) is an equilibrium solid 
solution of V in Fe-Co, in contrast to the supersaturated a, structure. 

The microstructures of Figs. 4 and 5 can be correlated with these 
diagrams (Fig. 6). At 950°C, the compositions are near the y phase 
boundary, hence the two-phase structure is primarily y with some ai 
(Figs. 4b and 5b).** The 3.0-percent V material is practically all y 


“The y phase formed by these anneals transforms to az on quenching identical to 
the transformation from the single-phase y region after the hot-rolling stage. There- 
fore, the two-phase microstructure is actually a1 + a2 at ambient temperature as 
photographed. However, the point of interest is the percentage of each phase formed 
at the elevated temperature which is identical to that existing at ambient for a 
quenched sample. Thus the designation y + a1 will be maintained for clarity in 
relating the microstructures to the phase diagrams. 

T The a1 etches preferentially to the a2(y) using the Nital etch. And, under the 
interference contrast conditions used in all cases, the a: phase (raised) appears 
bright on the top edge as if the light source were shining from the 12-0’clock posi- 
tion. Therefore, for example, in viewing Figs. 4b and 5b the primary or major phase 
should appear raised, whereas in viewing Fig. 5d the minor phase or particles should 
appear raised. 














/RIN/NI\IAX y 
fMLP 
A GOO 

















Fig. 6—Isotherms of the Fe-Co-rich portions of the Fe-Co-V ternary equilibrium diagram constructed from the vertical sections 
in the work of Késter and Schmid.! 
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(Fig. 5b) as this composition almost coincides with the phase boundary. 
At 900°C, the alloys are near the a; phase boundary, hence the two- 
phase structure is primarily a; with some y (Figs. 4d and 5d). In this 
instance the lower-vanadium (2.5-percent) alloy is apparently all a, 
(Fig. 4d) as this composition almost coincides with the a; phase bound- 
ary. For anneals at 925°C, these alloys should have nearly equal 
proportions of a; and y. This can be seen to be the case for the higher- 
vanadium alloy (Fig. 5c). However, very little y is present in the low- 
vanadium material (Fig. 4c). This may be an indication of a non- 
equilibrium character in these anneals. Apparently, the more highly 
strained ae lattice of the 3.0-percent V alloy aids sufficiently in the 
kinetics of nucleation of the y phase, such that this alloy approaches 
equilibrium more rapidly than the 2.5-percent V alloy. 
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925°C 
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Fig. 7—Quenched microstructures of 2.5% V-0.53-mm Remendur wire; etched ; 
DIC; 600X. First anneal indicates temperature of } hour anneal of 6.385-mm rod. 
Second anneal indicates temperature of 4 minute strand anneal of 0.53-mm wire. 
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A second influence of the higher vanadium content is apparent in 
comparing Figs. 4 and 5. A significant degree of grain size refinement 
occurs for the 3.0-percent V alloy. 

All six of the annealed microstructures were found to be ductile. 
They were drawn from 6.35-mm rod to 0.53-mm wire by successive 
20-percent reductions (23 passes) without any failures or evidence of 
cracking. This result clarifies several points with regard to ductility 
in Remendur. It has previously been shown that the air-cooled a, 
phase is brittle. English has shown that the ai phase orders (a) on 
‘other than very rapid rates of cooling? and alludes to the enhanced 
brittleness of this microconstituent.® Stoloff and Davies presented 
more complete evidence of the significant loss of ductility due to 
ordering.!! Chen’ and Davies!* have shown that some degree of ductility 
also can be created in the a2 constituent by the rapid cooling provided 
by a quench. However, Chen® demonstrated that the quenched single- 
phase a2 structure does not provide the degree of ductility of the 
quenched two-phase structure. The influence of the quench on the az 
constituent may also be related to the ordering phenomena although 
this has not yet been verified. Thus the two-phase structures developed 
by anneals between 900 and 950°C and followed by a quench are be- 
lieved to provide the maximum degree of ductility attainable in this 
alloy. Anneals at temperatures below 900°C but above the ordering 
temperature (720°C) followed by a quench may also produce a 
ductile structure. However, structures produced by lower-tempera- 
ture anneals for drawability may be undesirable for developing proper 
magnetic properties in the final reed. 


3.3 Annealed 0.53-mm (0.021-inch) Wire 


Remendur rods of both compositions and all three annealing tem- 
peratures were drawn directly to 0.53-mm wire (6 samples). Sections of 
each sample were strand-annealed at furnace set temperatures of 
900, 925, and 950°C for 30 seconds to produce a series of 18 samples 
with different combinations of composition, 6.35-mm anneal tempera- 
ture, and 0.53-mm anneal temperature. The actual peak temperatures 
as determined by drawing a thermocouple through the furnace were 
935°C, 913°C, and 890°C. The appropriate microstructures are shown 
in Figs. 7 and 8. The most significant points to be noted are as follows: 


(t) For the low-vanadium alloy (Fig. 7) the second anneal uniquely 
determines the percentages of a: and y. (Compare the identical micro- 
structures of each horizontal row.) This is due to the fact that, con- 
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Fig. 8—Quenched microstructures of 3.0% V—0.53-mm Remendur wire; etched; 
DIC; 600X. First anneal indicates temperature of } hour anneal of 6.35-mm rod. 
Second anneal indicates temperature of 4 minute strand anneal of 0.53-mm wire. 


trary to the case for the anneals of the 6.35-mm rod, the anneals in 
the 0.53-mm wire apparently produce a nearly equilibrium two-phase 
structure. This is likely a consequence of the aid to kinetics of nuclea- 
tion provided by the deformation of wire drawing. 

(iz) For the high-vanadium alloy (Fig. 8) the second anneal also 
primarily determines the percentages of a1 + y in qualitative agree- 
ment with the equilibrium diagrams (Fig. 6). The first anneals have 
little effect on phase percentages in the annealed wires indicating the 
0.53-mm, 30-second strand anneals may produce equilibrium structures. 

(iz) The same refinement in grain size of the high-V alloy compared 
to the low-V alloy previously noted is again present. An additional and 
substantial refinement occurs in the reduction of the wire from 6.35 
to 0.53 mm (compare Figs. 5 and 8). 
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3.4. 600°C-Annealed 0.538-mm Wire and 0.18-mm Ribbon 


Each of the 18 samples was cut into two lengths. One length was 
rolled to a thickness of 0.18 mm (0.007 inch) and the other length 
remained as 0.53-mm wire. These conditions simulate the paddle and 
shank, respectively, of a typical reed. All 36 samples were given the 
standard 600°C, two-hour anneal in dry H» which is used to produce 
a square magnetic hysteresis loop with a specified range of coercive 
force. 

The values of the coercive force developed for the various combina- 
tions of strand-annealing temperature and vanadium content in the 
round and flat sections are listed in Table II, and representative micro- 
structures for both round and flattened sections are shown in Figs. 9 
and 10. The details of the phase transformation which occurs during 
this 600°C anneal are presented in another paper! and will not be 
restated here. It is sufficient to note that the final microstructure 
developed in the round (undeformed) wire by this anneal is a strong 


TaBLE II—CoERcIVE ForcE IN OERSTEDS OF 0.53-mM ROUND WIRE 
AND 0.18-mMmM Fiat RIBBON IN THE 600°C—ANNEALED CONDITION 


(2.5 wt. % V) 
First Anneal 


950°C 925°C 900°C 
F R F R F R 
ta 
= 950°C 22.9 11.9 23.7 11.9 24.1 11.5 
~ 925°C 21.8 11.0 22.9 11.5 93.7 11.5 
i=] 
8 900°C 22.6 11.9 22.9 83 93.3 5.1 
D 
(3.0 wt. % V) 
First Anneal 
950°C 925°C 900°C 
F R F R F R 
a 
950°C 29.7 18.2 30.8 20.2 30.8 20.6 
* 925°C 30.8 23.7 32.0 24.9 32.0 24.6 
| 
8 900°C 30.5 17.0 30.8 17.8 31.3 21.0 
M 
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Fig. 9—Microstructures of 0.53-mm round wire and 0.18-mm flattened ribbon 
following 600°C, 2-hour anneal of 2.5% V alloy; etched; DIC; 1500X. Subcaptions 
specify temperatures of first and second anneals prior to 600°C anneal. 


function of the microstructure from the strand-annealed state. This 
may be seen by comparing Fig. 7 to Fig. 9 and Fig. 8 to Fig. 10 for the 
low- and high-vanadium alloys, respectively. For the low-V alloy strand 
annealed at 900°C (Fig. 9a) the large grain a matrix with particles 
of y at the grain boundaries produces a coercive force of only 5.1 
oersteds. The other structures show a coercive force of approximately 
11.0 oersteds with a maximum value of only 11.9 oersteds. This demon- 
strates the significant influence of vanadium content on coercivity and 
the need to keep the minimum permissible vanadium content above 
the 2.5-wt. percent level to achieve the required minimum value of 18 
oersteds on the round section of the reed for remreed applications. 
The high-V alloy shows a clear correlation between magnetic proper- 
ties after the 600°C anneal and the prior strand-annealed microstruc- 
ture. As shown in Table II a peak in coercivity of the round wire is 
attained for the 925°C strand anneal of 0.53-mm wire. This is the 
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Fig. 10—Microstructures of 0.53-mm round wire and 0.18-mm flattened ribbon 
following 600°C, 2-hour anneal of 3.0% V alloy; etched; DIC; 1500X. Subcaptions 
specify temperatures of first and second anneals prior to 600°C anneal. 


state of a nearly equal distribution of the two-phase a1 + az structure 
(Fig. 8) which produces a more uniform distribution of a, (ordered a1) 
and y (Fig. 10b) during the subsequent 600°C anneal. Coercivity de- 
creases for strand anneals at higher or lower temperatures where the 
a1 OF a2 phase predominates. 

The flattening operation eliminates most of the microstructural 
influence on final magnetic properties as shown in Table II. For the 
high-V alloy, all values are in the range of 30 to 32 oersteds and for the 
low-V alloy in the range of 22 to 24 oersteds. All microstructures are 
similar to that shown in Figs. 9d and 10d which consists of the very 
fine distribution of y in an a matrix. 


IV. CONCLUSIONS 


The primary conclusion from this study is that the microstructure 
of 2-3-wt. percent vanadium Remendur is extremely sensitive to 


PROCESSING OF REMENDUR 13389 


annealing temperature and vanadium content. Magnetic characteristics 
of Remendur wire and the remreed contact and the mechanical proper- 
ties, particularly drawability, of the material are dependent on the 
precise microstructure developed during processing. More specific 
conclusions are listed as follows. 


(t) Anneals in the range of 900 to 950°C produce a two-phase 
ai (BCC) + vy (FCC) structure which changes to a two-phase a; (BCC) 
+az (supersaturated BCC) structure on rapid cooling. The percentages 
of a; and y formed are a strong function of temperature with increasing 
temperatures yielding greater percentages of y (a2). 

(2t) The one-half-hour anneal of low-vanadium 6.35-mm rod does 
not produce the equilibrium two-phase structure defined by the phase 
diagram. Increased vanadium content assists in producing a nearly 
equilibrium structure as does the deformation of wire drawing for the 
anneal of 0.538-mm high- and low-vanadium wire. 

(i727) The percentages of a; and y (a2) formed in the 0.53-mm wire 
are primarily a function of the temperature of the anneal. The tem- 
perature of the 6.35-mm rod anneal has little influence on determining 
these final percentages. 

(tv) The all ae structure created by slow to moderate rates of cool- 
ing from the y phase region is brittle and cannot be cold drawn. How- 
ever, a two-phase mixture of ai + a2, developed by quenching from 
anneals in the range of 900 to 950°C (the two-phase a + y region), 
is sufficiently ductile to be drawn from 6.35-mm rod directly to 0.53- 
mm wire. If the rate of cooling is not sufficiently rapid from the 900 to 
950°C range, the a; structure orders? to form a; which has also been 
shown to be brittle and, therefore, nondrawable to wire." 

(v) Increased vanadium content and the deformation of cold drawing 
both are influential in refining the grain size of the two-phase Remendur 
structure. The refined grain size and higher vanadium content in- 
fluence the magnetic properties of 0.53-mm Remendur wire which has 
been annealed at 600°C for 2 hours by increasing the coercivity. 
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Optimum Mean-Square Decision 


Feedback Equalization 


By J. SALZ 
(Manuscript received April 20, 1973) 


In this work we report new results relating to decision feedback equaliza- 
tion. The equalizer and the transmitting filter are optimized in a PAM 
data communication system operating over a linear noisy channel. We 
use @ mean-square error criterion and impose an average power con- 
straint at the transmitter. Assuming correct past decisions, an explicit 
formula for the minimum attainable mean-square error is given. The 
possible advantages of signaling faster than the Nyquist rate while de- 
creasing the number of levels to maintain the same information rate are 
investigated. It is shown that, in all cases of practical interest, signaling 
faster than the Nyquist rate, while keeping fixed the information rate, 
increases the mean-square error. Finally, to illustrate the use of the 
results, application is made to a cable channel where the loss in dB varies 
as the square root of frequency. Various asymptotic formulas and curves 
are provided to exhibit the relationships between the quantities of interest. 


I. INTRODUCTION 


A great deal of research, particularly in the past decade, has been 
expended on the problem of linear equalization. This has yielded a 
considerable body of theory and technology making possible the design 
of apparatus for successfully combating intersymbol interference in 
PAM data transmission systems operating over noisy linear channels 
where delay distortion predominates. Since linear equalizers must com- 
pensate for the channel characteristics in the presence of noise, they 
cannot be expected to perform well over severely frequency-attenuating 
channels or channels possessing nulls in the amplitude characteristic. 

Interest in the high data rates over voiceband and cable channels 
inevitably leads to the search for more effective equalization methods. 
Faster pulse rates place signal energy well within the badly attenuated 
portion of the transmission spectrum, resulting in severe intersymbol 
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interference correctable by linear methods only at the expense of a 
significant enhancement of the noise. 

A “bootstrap”? technique, commonly referred to as ‘decision feed- 
back,”” when combined with linear equalization can yield significant 
performance improvement.!? In this method the samples of the pulse 
tails (postcursors) interfering with subsequent or future data symbols 
are subtracted without incurring a significant noise penalty. The effect 
of pulse tails (precursors) which occur prior to detection and inter- 
fere with past symbols is minimized by a conventional linear equalizer. 

Much has been written about this subject. In a fundamental paper 
where an excellent bibliography can be found, Robert Price? demon- 
strated quantitatively the merits of decision feedback equalization in 
certain applications. 

In this work we jointly optimize the receiving and transmitting 
filters in a PAM data transmission system employing decision feed- 
back. The chief difference between our work and Price’s is in the choice 
of performance criterion. We minimize mean-square error while 
Price maximized the signal-to-noise ratio under the constraint that the 
overall intersymbol interference be zero. Our criterion is not as strin- 
gent as Price’s and allows trade-offs between added noise and inter- 
symbol interference. Monsen‘ also investigated some aspects of our 
problem but did not arrive at a complete solution. Our chief contribu- 
tion is an explicit formula for the minimum mean-square error (MSE). 
The simplicity of the formula makes possible detailed investigation of 
optimized system performance. 

In Section II the model is stated and the problem is formulated. 
In Section III the receiving filter is optimized and in Section IV the 
transmitter filter is optimized. In Section V we examine the problem 
of signaling faster than the Nyquist rate and finally in Section VI 
we use our results to investigate in detail the performance of a data 
system operating over a cable channel. 


Il. THE MODEL AND PROBLEM FORMULATION 


The system model under investigation is depicted in Fig. 1. The data 
signal denoted by D(t) is passed through the transmitting filter having 
an impulse response s(t) and giving rise to an average transmitted 
power P. The data symbols {a,}%. are independently picked at the 
rate 1/T and take on values with equal probability from the set 
{t1, +3+45---+(L—1)} where L is an even integer. The 
resulting signal is admitted to a linear channel characterized by an 
impulse response h(é). The received signal plus noise is processed by 
the equalizer which is comprised of a linear filter having impulse 
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Fig. 1—Block diagram of the model. 


response w(t), a sampler, a decision rule, and a feedback digital filter 
characterized by the infinite set of real numbers {b,};°. The added 
noise n(¢) is a zero-mean white random process with double-sided 
spectral density No/2. The output data symbols are denoted by 4n. 

The general problem we would like to solve is the minimization of 


MSE = E{v, =x Gx}? 


with respect to the set of square integrable functions {s(é), w(t)} and 
the infinite sequence of numbers {b,}. This is to be carried out when 
the channel impulse response h(é), transmitted power P, and a decision 
rule are given. The symbol #(-) denotes expectation with respect to 
all the random variables. 

The nonlinear relation between the estimated symbols {dé,} and the 
input symbols {a,} makes this problem mathematically intractable. 
However, by assuming that past decisions have been correct, we can 
begin to approach the problem. The resulting MSE must then be 
interpreted as a lower bound on the true MSE. Alternatively, if no 
errors have occurred in the past, the MSE under this assumption 
provides an indication of the noise immunity of the system (including 
residual intersymbol interference). 

Let r(t) = s(t)*h(t)*w(t) denote the overall impulse response, 
where * denotes convolution. Under the assumption of correct past 
decisions, the received sample taken at time ¢ = kT is 


ioe} 


= DY Prdz—n — x brdz—n + n(t) *w(t) | mar 


n=—0o 


and the mean-square error is then by definition 


—Il rove) 
n % n=1 


Sb nen heehee Wael 
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A straightforward calculation gives 


= 


—1 00 
MSE = o DY r24+ 62 (r) — 1)? +62 > (12 — bn)? + 0’, 
re) n=] 


where 
i—1 
3 


oa a fe w(B)db. 


It can be immediately concluded that the mean-square error is 
minimized by setting bn = rn, n = 1, 2, ---, which eliminates the 
feedback coefficients {b,} from further consideration. 

The problem we now confront is the dual minimization of 


= Kia}? = 





and 


MSE[s(Z), w(t)] = o| D2 tlt =| 


with respect to s(é) and w(t) when a constraint is imposed on the 
average transmitted power. 

The above expression indicates that, under the assumption of 
perfect past decisions, the mean-square error is minimized by minimiz- 
ing both the pulse precursors in the overall impulse response and the 
output noise power, while keeping 7» close to unity. 

We give a precise formulation and solution to this problem in 
Section ITI. 


III. RECEIVER OPTIMIZATION 
Writing the MSE in detail we obtain 


Meni t & | [weer — oar| 


a 





91 an =Dde +33 c Gdn. A) 


—* 


where 
p(t) = s(t) *h(t). 


Keeping s(é) fixed, and using a standard calculus-of-variation ap- 
proach, results in an integral equation for w(t) 


0 0 
r(-) = Nw + [wi (par —np@r—o) dr, @ 
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where 
202 


No 


If in eq. (2) we set 


ie ie w(r)p(nT — z)dr, 


we see that the optimum solution must have a representation in the 
form 


0 
w(t) = ms gnp(nT — t), (3) 
where 
1 
go = ip Ch) 
and 
Un 
n= yy? 29 n= s = 1 
g Mm, 


In (8) is revealed the structure of the optimum receiving filter. It is 
composed of a matched filter having impulse response p(— t) followed by 
a one-sided (anticausal) tapped delay line with weights equal to gn. 

Linear equations involving the set {U,} can be obtained by first 
multiplying both sides of (2) by p(kT — t), k S$ 0, then integrating 
from minus to plus infinity. The resulting linear system of equations is 


0 
R, = N\Uz + Kia 4»; k = 0, =] sx, (4) 


where 


Ree [ p(— t)p(kT — t)dt. 


The system of eqs. (4) can be solved by standard Wiener-Hopf tech- 
niques and the details are given in Appendix A. The solution in terms 
of the discrete Fouricr transform of the sequence {U,}°. 1s 


0 
Ss U,ei”?® 


n= oe 


U(6) 


N 
1 IE Or is 
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where 
M(6) = M+(0)M-(¢) = R@)+N,= x M,ein?, 
M+ (@) = > yne'”9, 
n=0 
and 


M-(6) = M+(— 6). 


This is standard procedure making use of the well-known factoriza- 
tion property of covariance functions. Methods for obtaining the 
sequence {y,}¢° from the given sequence {R,}2. are well documented 
in the literature. One method is summarized in Appendix A. 

Having specified the optimum receiving filter, we now obtain a 
formula for the minimized mean-square error. The availability of this 
simple formula will allow further optimization of the transmitting 
filter. 

Let wo(t) be the impulse response of the optimum receiving filter. 
[This function solves the integral equation (2). ] Substitute wo(t) 
into (2), multiply both sides by wo(t), and integrate from minus in- 
finity to plus infinity to obtain 


[i (— dwo(pat = No [wear 
+ & (fer - Dwo(dat)~ (6) 


n=—-2 


Putting this into (1) with w(t) replaced by wo(t) we get a formula for 
the optimized MSE 


MSE[wo(t)] = o&(1 — Uo) = Nooago. (7) 


This result was obtained by Monsen‘ but unfortunately he did not 
go any further. As it turns out, a much richer formula than (7) can be 
obtained since Up can be expressed directly in terms of the spectrum 
of the channel characteristics in cascade with the transmitting filter. 
To carry this further, observe from (5) that 


Uy = de term in U(@) 
No 


es fee 
v5 


(8) 





and consequently 
. Ni ee 


v5 


MSE = «2 (9) 
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As it turns out, yo is functionally related to M(@) in a rather simple 
manner. This relationship can be found in the literature but the deriva- 
tion is short and so we briefly outline the approach. 

Under very mild conditions on M(@) (see Doob,® pp. 159-161) 


M(8) = : (10) 








DD ynei"? 
n=0 
where yo is real and positive. Since (@) > 0, consider 


" In M(6)dt = ic iA E iS vaein [do 
—7 n=1 


i a In | vo + Xe la. (11) 


When the In’s on the r.h.s. of (11) are expanded in a power series 
and the integrations are carried out (recognizing that all integrals 
involving powers of exp {in@}, n ¥ 0, vanish) we get 


= exp | 2 [7 InlR@ + Nol, (12) 
where 
R() = X Rret, 
- : dw 
= 2ptwnT 
Ry =f |Plw) |r? 5, 
and 


PG) = [ : s(t) #h(t)e**dt. 


After minor algebraic manipulations and changes of variables we 
obtain by substituting (12) into (9) 


MSE = ¢? Pf tary 1d 13) 
g = otexp |—5- [int (0) + 1]do!, —— ( 
where 
] Qian \ |? 
Yo) = rp P(o- 7) (14) 








This formula, as far as can be determined, is new and its simple 
form will enable us in Section IV to carry out an additional optimiza- 
tion with respect to the transmitting filter. 

It is instructive at this point to compare this formula with the one 
obtained for a linear equalizer without decision feedback. Berger and 
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Tufts® have found such a formula and from their paper we have that 


r/T 
(MSE)tinear = 02 5 [" [¥(@) + Udo 
W J—a/T 


= aei) (8) 


[T priv 
. me Y aa . ad . 
() cis. fC: Me 


where 


In terms of the same notation, (13) can be put into the form 
MSE = exp{—(Un[Y() + 1])} (16) 


from which we get immediately that 


MSE < (erwin) = (aoe) = (MSE) inear- (17) 


As expected, the mean-square error with decision feedback is always 
smaller than the MSE of a linear equalizer. Comparing (15) with (16) 


1.00 


AMPLITUDE SQUARED 
° 
a 
{fo} 


0 1000 2000 3000 4000 5000 
FREQUENCY IN HERTZ 


Fig. 2—Amplitude-squared characteristic of typical voiceband channel. 


DECISION FEEDBACK EQUALIZATION 1349 


—10 LOG N’9=~40 dB 


LINEAR EQUALIZER 


MSE IN DECIBELS 


FEEDBACK EQUALIZER 





0 2000 4000 6000 8000 10,000 
RATE IN BITS/SECOND 


Fig. 3—MSE in dB vs binary data rate for channel shown in Fig. 2 without de 
transmission. 


shows that both equalization methods yield the same MSE if and only 
if Y(w) is a constant, i.e., there is no intersymbol interference. 

Prior to optimizing the transmitting filter we wish to illustrate the 
behavior of (15) and (16) for a typical voiceband channel as the signal- 
ing rate 1/T is allowed to increase. An amplitude-squared character- 
istic for a typical voiceband telephone channel is shown in Fig. 2 
(the dashed line with zero transmission at zero frequency is typical). 
Figure 3 shows the resulting MSE vs pulse rate for both a linear equal- 
izer and a decision feedback equalizer when o2 = 1 (binary data) and 
— 10 log Ny 40 dB. The calculations were done numerically by 
using (15) and (16). We note that the performance of the linear 
equalizer deteriorates rapidly when the rate is greater than = 3000 
bits/second while the decision feedback equalizer deteriorates grace- 
fully. The reason the linear equalizer shows such a poor performance 
is that it must compensate for the missing energy around dc. In 
practice, of course, modulation is used to place the data energy at a 
more suitable location in the passband spectrum to avoid this severe 
null. To make the comparison fair, we artificially extended the channel 
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—10 LOG No = 40dB 


LINEAR EQUALIZER 


MSE IN DECIBELS 


FEEDBACK EQUALIZER 





0 2000 4000 6000 8000 10,000 
RATE IN BITS/SECOND 


Fig. 4—MSE in dB vs binary data rate for channel shown in Fig. 2 with de 
transmission. 


characteristic from about 300 Hz to 0 Hz to a constant transmission. 
This is indicated in Fig. 2 by the dashed line parallel to the frequency 
axis. Figure 4 shows the comparisons for this atypical channel. As 
expected, the linear equalizer has a sharp threshold at approximately 
the Nyquist rate but, as before, the decision feedback equalizer 
deteriorates much more gracefully. 


IV. TRANSMITTER OPTIMIZATION 


The problem we address here is the optimization of (13) with 
respect to the transmitting filter characteristics subject to an average 
power constraint. 

Let S(w) and H(w) be the Fourier transforms of s(t) and h(t) re- 
spectively. The average power at the output of s(f) is 


2 ) 2 oe) 
pS Z| (dt = | S2(w) do 


Zea ip S ged een 
= Foe [= 8 (e-ae)) a a) 


where by S?(w) we mean |S(w) |?. 


i 
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The problem at hand is to maximize the functional 
r/T 
I= i In[K ¥ S2(w)H2(w) + 1]de 
—xr/T n 
w/T 
+A LU Si(w)jdw (19) 
/T n 


with respect to the infinite set of functions {S2(w), n = all integers}. 
Where 
S2(e) = S(c — a ) 


and \ is a Lagrange multiplier to be determined from the constraint 
on the average transmitted power. Observe that these functions 
are independent over the range — 7/7 Sw S7/T and therefore 
consider the variation of J with respect to, say, S3(w) 


2 


. K=1/NiT 





ees r/T KH?(w) 2 ° . - 
6,1 [ - KEBAB) i BSH) + MSH(w) | deo. (20) 


When S37(w) # 0, setting (20) to zero implies 


KH3(o) : 
Kx S@)Hi~) tot *=° = 
for all wC€ [—7a/T, x/T]. Now, assume that H?(w) > H2,(w) for 
|n|<|m| and » € [—2/T, 7/T ]. When this condition is satisfied it 
is not possible to solve the system of equations given in (21) unless 


S2(w) = 0 for alln ¥ j 
in which case we get 
KH3(w) = — _KS3(w)Hj(w) + 1]. (22) 


Substituting this into (20) indicates that the largest value is obtained 
when S3?(w) = S2(w) = S?(w) and furthermore \ must be negative. 

So far we can conclude that for channels possessing monotonically 
decreasing amplitude characteristics, i.e., H?,(w) > H2(w), — «/T 
Sw S7/T, when |n|>|m|, the optimum transmitting filter cuts 
off at the Nyquist frequency 7/7’. The optimum system allows no 
transmission outside the band |w| > 7/7. The restrictions imposed on 
the channels are mild and are expected to be satisfied in most situa- 
tions of interest. However, removing these restrictions makes the 
problem slightly more complicated, and it is left up to the reader to 
reason how it can be solved. 
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Making use of this partial solution permits writing the mean- 
square error in the simplified form 


E } Oe 


Sin (“> ) a [om [KH?(w)S2(@) + 1]do (28) 


and the functional now to be further maximized with respect to the 
inband structure of S?(w) reduces to 


ise) = [ “In [KS2()H2() + Udo + a  Sto)de, (24) 


As can be seen from (20) and (22), eq. (24) is maximized when 


Kiley 253 0] 


\KE2(a) (25) 


S?(w) = max | 
where \ = — Xo. 
To determine the Lagrange multiplier \9, two cases must be dis- 
tinguished. 
Case 1: KH?(w) — Xo > O, for all |w| S +/T 
In this case, we get 


1+ KH?S? =yH?, w= K/r, 


which when substituted into (23) results in an explicit expression for 
the optimum MSE 


a/T 
—In (=) =Inn+— [Un Pode. (26) 
0 


a 





The factor u is determined from the average power constraint in the 
following manner. Use (25) and (18) to write 


_o% f77[u_ 1 
Por |, EF Kw | 
4 
= 77 (# — A), (27) 


where 


T pur 4 
A==f Po %* 


Substituting the parameters K = 1/TN 5 and Ng = No/2c2 into (27) 
gives explicitly 
w=pt+A, (28) 
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where 
<2 re _ Average transmitted signal power 
p ( No 7) ~ Average noise power in the Nyquist band 
yagi 


Thus (26) and (28) provide a complete solution for this case. 


Case 2: There exists a set of w for which KH?(w) — \»o S 0 


The optimization procedure in this case involves the standard water- 
pouring argument. To illustrate the nature of the solution we take 
the situation where H?(w) is strictly monotonically decreasing in the 
Nyquist band. This implies that there exists only one frequency wo 
for which KH?(w) = 9 and consequently we get 


1 
"FP (ao) 
This gives one relation between the unknowns y and wy and another is 


obtained from the power constraint. Since the optimum filter char- 
acteristic is zero when w > wo, the signal-to-noise ratio 1s 


(29) 





T so 1 1 
ee l Fay ~ Fey |” (30) 
and an explicit formula for the mean-square is 
= In | Mee | po r In FH? (w)dw = Fw In H? (a0). (31) 
a 0 


We now briefly summarize how these optimized formulas are to 
be used: 


(<) For a given transmitted average signal-to-noise ratio p, solve 
eq. (30) for wo. 
(22) If wo < 7/T, use formula (31) to compute MSE. 
(272) If wo 2 +/T, use formula (26) to compute MSE. 


In Section VI we shall illustrate numerically the use of these 
formulas. 


V. SIGNALING FASTER THAN THE NYQUIST RATE 


Here we examine the behavior of the optimized mean-square error 
when the frequency support of an ideal unity-gain channel is smaller 
than the Nyquist rate $7. After deriving the optimized MSE for this 
situation, the possibility of further optimization relative to the signal- 
ing rate when the information rate per unit bandwidth is held fixed 
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will be investigated. This has been an open question thus far and the 
issue is whether increasing the signaling rate beyond the Nyquist rate 
while decreasing the number of levels to maintain a fixed information 
rate is ever beneficial. 

Consider a channel having the characteristic 


0, we By 


H?(w) = 1 w ‘= EB,’ 


(32) 
where the sets £; and HE, form a partition on the frequency interval 
E = {w:0 Sw =7/T}. By frequency support we mean the measure 
of the set H2 denoted by m(F2). 

For this channel it is easy to calculate explicitly the mean-square 
error. Observe from eq. (22) that for a piecewise constant channel the 
optimum transmitting filter is a constant when w C E>, and zero other- 
wise. The minimum MSE is then calculated from (23) 

—In | | ae ie In [KS?H?(w) + 1]de 
o 2a J_x/T 


a 


cod 2 
- i, ap, BLKS* + 1d 


= F m(Bs) In [KS? + 1), (33) 
where S? is a constant to be determined from the power constraint 
2 2 
= Sa 2 = Ga 2 
7 oe dw ps m(E2). (34) 


From this equation and the definition of K = 1/(N oT) it can be 
checked that 
p= KS = Average signal power 


Average noise power in a band = m(E2)- 
Substituting this into (33) gives the simple desired formula 


MSE = of(1 + p)7,; (35) 
where 


Tv 


_ Channel bandwidth 
~ 2X Signaling rate 


For fixed p, MSE — o? when a=0 and, as expected, MSE — o2(1+ p)™ 
when a = 1. It is curious that as long asa + 0, MSE—0 as p— ~. 


1 


IIA 
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When the set EH is the interval J = {w:0 Sw Sa/To < 1/T}, it 
is referred to as being less than the Nyquist band. Since there is no 
mathematical reason for making this distinction we shall refer to 
“signaling faster than the Nyquist rate’ whenever a < 1. 

We now investigate whether it is ever advantageous to signal 
faster than the Nyquist rate. Clearly, for fixed o2 = (L? — 1)/3, or a 
fixed number of levels, (35) shows that MSE degrades rapidly with 
decreasing a. An interesting question, first raised by R. W. Lucky,’ is 
the possibility of trading L with a to further minimize the mean- 
square error. This is the problem we address. 

Let the source information rate be R = logeL/T bits/second. 
The available bandwidth is equal to 1/(27)m(E2) cycles/second. Thus 
the normalized information rate is 


Pe 


= -,—, 5 _logel 
=~ m(B») ame) | 
T T 
2 bits 
Writing (35) in terms of this quantity gives 
226 — J 1 
MSE (a) = oa aa (1 + p= (37) 


Letting C = loge(1 + p) be the ultimate attainable rate according to 
Shannon’s theory, eq. (37) can be put into the form 


MSE(a) = (2-2(6-8) — Q-acy2, (38) 


Note that, since L = 2%8/? and the minimum allowable L is equal to 2, 
the parameter a must be in the range (2/6, 1). 

The problem initially posed can now be stated as follows. Find a 
set of a’s, 2/6 = a = 1, which minimize eq. (38). We begin by setting 
the derivative of (38) to zero, 


a = — (C — 8)2-#C-8) 4 CQ-aC = 0, 
a 


from which we find a unique stationary point 


2 Ci 
a= Zlogs| os | (39) 


Since MSE(O) = 0 and MSE(@) > 0, the value of a found above 
must be a point where MSE attains a maximum. If this maximum is in 
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the range (0 =a S 2/8), then the minimum value of MSE is attained 
at the boundary (a = 1). The condition for this to be the case is 
determined from (89). 


C PL 
logs loa | = 1. (40) 


From this we deduce that as long as 6 < 2C, a = 1 minimizes MSE. 
For this region of & and C there is no advantage gained by signaling 
faster than the Nyquist rate. Since C is channel capacity, & = 2C 
seems to be a critical rate. If the rate is greater than this critical rate, 
the maximum point lies within the allowable range of a(2/6, 1) thus 
raising the possibility that eithera = 2/& ora = 1is a minimum point. 
Suppose a = 2/6 is a minimum point. Equation (38) gives for this 
case 


) = 2-2C/8 > 9-813 0.157. (41) 


MSE (2 ~4 


Thus we have found a region for which signaling faster than the Ny- 
quist rate appears to be beneficial. But at this high level of MSE, we 
are no longer justified in assuming that the feedback decisions are 
correct most of the time. In fact, what is more likely to happen is 
that errors begin to occur resulting in a larger value of MSE than 
predicted by the error-free model. We are therefore led to the con- 
clusion that a minimum point other than at a = 1 will render an MSE 
to be outside the range of practical utility. To emphasize this point 
further substitute (41) into (38) to get 


MSE(« = 1) = [(MSE(2/8)J¢? >. 


(42) 


Thus a = 1/6 is a minimum point whenever 





(MSE (2/68) Je” (= - ~) > MSE(2/68) 


or 





3 
MSE(2/8) > (5 = ; Zo 

As an example of the use of these inequalities, suppose that & = 3 
and C = 3.5; therefore, a = 3 is the minimum point and the achievable 
MSE = 0.198 which can be obtained with a binary system (L = 2). 
On the other hand, M(a = 1) = 273-5(7/3) = 0.206 which can be 
achieved with (L = 2!.5 = 2.8)!! 

It is interesting to see what MSE can be achieved when only a 
linear equalizer is used. Using formula (15) and following the same 
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reasoning as before yields an expression for the optimized mean- 
square error 


(i eje tl 
. — 7 ip. ce i A 
(MSE) inear Ca p l 1 . 


In this case, unlike in the decision feedback case [eq. (35) ], as p—> ~, 
MSE — 1 — a, when a > 0. Thus the mean-square error cannot be 
made vanishingly small as the signal-to-noise ratio increases without 
bound. 

Expressing (43) in terms of the normalized rate & we get 


2 a, pia a 
(MSE) iinear — — = 5 | a 2 (44) 
Again we seek to minimize (44) with respect to a in the range 
(2/6, 1). It can be checked that (44) has at most one stationary point 
in the range 0 =a = 1.Sincea = 0isa minimum point, the minimum 
in the range (2/6, 1) must lie on the boundary. Thus the condition for 
achieving a smaller MSE when signaling faster than the Nyquist rate 
(a < 1) is 


(43) 


Pa | 


pl — 2/6) + 1S 3 





or 
2§ — I 1 

3 1-2/8 
Is it possible to find an 6 =C, 6 > 2 which satisfies the inequality 
(45)? A straightforward analysis reveals that the answer is negative. 
In other words, the Nyquist rate is optimum provided the information 
rate is less than channel capacity. 


p=2°-185 (45) 





VI. APPLICATION TO A CABLE CHANNEL! 


This section will illustrate the use of the formulas developed in 
previous sections in a particular application. For this purpose we 
choose a cable channel having frequency characteristic 


H(f) = exp {V— 2iaf}. (46) 


We shall develop in detail the applicable formulas, provide asymptotic 
behaviors, and exhibit numerically the relevant parameter trade-offs. 
For comparison purposes, the applicable formulas for the optimum 
linear equalizer will also be developed. 

tI am indebted to Dr. Robert Price for calling my attention to Ref. 8 where 


related work is reported. The paper is in Japanese; however, Dr. Price has an English 
translation. 
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We begin by first considering a suboptimum system where the 
transmitting filter is flat across the Nyquist band and zero outside. 
For this case, the minimum attainable mean-square error is [eq. (23) ] 


nw/T 
M, = exp |- f / In [S*Ke-VWeel* +. 1 |dw}, (47) 
0 
where 
M,; = MSE/c?, | 1 (w) |? = e~Veeeln, and K = 2¢2/TNo 


are to be determined from the average power constraint. Since the 
transmitted power P = (¢2/T)S?, we find that S?’K = 2PT/N o = p, 
the transmitted signal-to-noise ratio, where the noise is measured in 
a band = 1/7. Substituting these constants into (47) and making 
some changes of variables result in 


4 
Ma = exp |- 2f In [pe-V484 + 1 ]dy;, (48) 
0 


where 8 = a/T'. The parameter ¥28 is seen to be proportional to the 
loss of the cable in dB at the Nyquist frequency $7. 

We are interested in the behavior of (48) when p and 8 are varied. 
While it is not possible to express this integral in a closed form, it is 
possible to obtain a rapidly converging series in the two parameters of 
interest from which asymptotic behaviors can be deduced. Appendix B 
shows the details of the development. Different power series apply in 
different regions. The first series applies when In p < V¥28 and the 
second In p = v2@. The results are as follows 


1: V¥28 =Inp>0 


In p)® + 7? In 2 = e-V28]n 
—InM, = Serer ne _ ne pyres Loe 


Le Ea ERs 
+5 o (-S| 5 - oo]. (49) 


2: In p = V28>0 


0 26 Jn 
— In Ma = Inp ~ 3 198 + Res (— pe a(S | 


Hg coms[s-(S)] © 


It can be checked that (49) equals (50) when In p = V26. The first 
asymptotic behavior is deduced from (50) when p— © and 26 is 
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held fixed. For this case we get 


428 
os, (51) 


ay 





Another asymptotic behavior is deduced from (49) as B > ~ while 
p is held fixed. In this case 


Ma ~ e918, (52) 

where 

_ (Un p)? + 7? ln p ay Se er ee 

g(p) = 6 ge 2 ae bert 
In order to make possible performance comparisons, we develop 
similar formulas and asymptotes for a system employing only a linear 
equalizer. The minimum mean-square error applicable in this situa- 
tion is obtained from eq. (15). After substituting the cable char- 


acteristic we obtain 


5 
My = 2 f° [eV + 1p, (53) 
0 
where 
Mr = See 


Here, as in the decision feedback case, rapidly converging series can 
be developed from which asymptotic formulas are deduced. The 
detailed calculations are also given in Appendix B. The desired results 
are 


pais ages 


+5h(- tf + (pe¥8)*]. (54) 


=e 


2:In p > V¥28>0 


Mz = [2m ji+ | 


+35 (= yen | 5 (=)']. (55) 


p” 


The asymptotic formulas are readily deduced from (54) and (55). 
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When p — © and @ is fixed we get 


Mz ~ Q(8) -, (56) 


Q(B) = ev 2 -3+ > |. 


On the other hand, when 8 — ~ and p is kept fixed, 


where 


3(In p)? + 2?  D 
Mpw~] = —, 57 
i 6p Ts (57) 
where 
a 1 
— _ nt+1 — 
D= Z(- DSS. 


At this point it is possible to make some judicious performance 
comparisons between the two schemes. The first observation is that, 
in order to get a small mean-square error, p must be large and greater 


than e¥?6, In this case, (51) and (56) apply. On the other hand, when 


428 >In p, and 8 ~, performance deteriorates rapidly as can be 
seen from (52) and (57). Suppose now that a large signal-to-noise 
ratio is available and we wish to obtain the same mean-square error 
in both schemes. How do the signaling rates compare? 

Equating (51) and (56) shows that Ba/8z ~ 9/4 when these quanti- 
ties are large. In other words, asymptotically, the signaling speed of 
the cable may be increased by more than a factor of two with the 
use of decision feedback. Clearly when @ is small no significant ad- 
vantage can be obtained from using decision feedback equalization. 

To exhibit these phenomena further, we have used numerical 
integration to evaluate (49) and (53) and checked the accuracy by 
summing terms in the various power series. The results of these 
calculations are exhibited graphically in Figs. 5 through 9. A striking 
feature in all these curves is the manner MSE degrades as @ increases. 
The linear MSE exhibits a sharp threshold while the MSE for the 
decision feedback equalizer degrades much more gracefully. 

Next we wish to examine the possible payoffs when the inband 
characteristics of the transmitter filter are optimized. To do this 
explicitly, we follow the procedure outlined in Section IV. Equation 
(30) must first be evaluated for the cable characteristics. (We omit 
all straightforward integrations and algebraic manipulations.) Re- 
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_- V2Bo= 4.99 
a 


—10 LOG MSE 


_ DECISION 
// FEEDBACK 


V 2B = 6.35- — 


“ 
LINEAR — ~~ 





Vv 2p 
Fig. 5—MSE in dB vs ¥28 for p = 20 dB. 


writing eq. (30), 


= ery = Lawcea — ars | 


2 
_ Bol vee, _ F (280) | 
BLS Bo Pe 
where 
Bo = a/T o, oo = m/To, B= a/T, 
FH? (cw) = e~Vooal=, 
and 


F(x) = e¥#(vx — 1)+1 ~ 2/2 when z is small. 
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_U— V2Bo= 9.421 


—10 LOG MSE 
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V 2B, = 11.246-— 


a 
LINEAR ~ 





Fig. 6—MSE in dB vs V2¢ for p = 40 dB. 


The explicit evaluation of MSE is as follows: For a given p, 8 solve 
(58) for Bo. If Bo > B, calculate MSE from eq. (26), 





n/{[T 
—In || =InMs= Inu += [™ In H%(e)de. 
a 0 
An explicit evaluation gives 
328 
Me ee: > p: 59 
par =e 


On the other hand, if (58) yields a Bo < 8, use formula (31) to compute 
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_ —V 2Bp = 13.958 
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V/ 289 = 16.108 ——~ 
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V2B 
Fig. 7—MSE in dB vs ¥28 for p = 60 dB. 


MSE. The explicit evaluation for this case gives 


Mg = 67366018) v280, Bo > B, (60) 
1 4(Bo/B) 
p 5 + F(260)/Bo 
0 


where e¥?60 was obtained from (58). It can be checked that when 
Bo = B in (58), eq. (59) equals (60) as it must. 

Let us now pause and examine what these optimized results are 
telling us. Suppose @ is fixed in (58) and p is allowed to increase. 
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Vv28 
Fig. 8—MSE in dB vs ¥2¢ for p = 80 dB. 


Eventually a 8» will be found which satisfies (58) and which ultimately 
will be greater than 6. The physical implication of finding a 8» which 
is less than 8 is that the transmitting filter cuts off before the Nyquist 
frequency 47. This will occur only when p is relatively small and thus 
results in a poor MSE. Practically, the region of interest is when p 
is large such that By = 8, in which case the filter cuts off at the Nyquist 
frequency. In this region (59) applies and, upon comparing (59) with 
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p = 100dB 


frp = 23.1 


——V 2 Bo = 23.1 B = 
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—10 LOG MSE 
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LINEAR —~ 
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V2 
Fig. 9—MSE in dB vs V2¢ for p = 100 dB. 


the asymptotic suboptimized result, (51) shows that 
etv28 et V28 


p +=) p 


A 


Since ming [F(6)/8] = 1, the optimized result appears to be 
asymptotically equal to the suboptimized result. In other words, in 
the region where In p > ¥28 and p > ~ no benefits are obtained from 
inband optimization. This is also evident from eq. (25) since when In p 
is large relative to V¥28 the structure of the optimum transmitting 
filter is a constant. The situation where Bo < # is slightly more com- 
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plicated to compare. Here inband optimization should perhaps be 
beneficial. However, comparisons in this case between the optimized 
MSE and the suboptimized ones must be made on the basis of the 
same transmitted power rather than signal-to-noise ratio because the 
systems operate over different bandwidths. 

Again for comparison purposes, we summarize the formulas that 
apply when the inband characteristics of the transmitting filter in a 
system using only linear equalization are optimized. Berger and Tufts® 
carried out such an optimization and the procedure is similar to the 
one carried out in Section V. Adopting our notation and the same 
definition of parameters as above we can obtain explicitly the follow- 
ing formulas applicable for a linear system. 

Choose a p and a @ and solve for 8» in the equation below: 


_ Bol 4 vin — 1 
p = 5) 5 F(G0/2)e¥# — 3 F (280) |. (61) 
If Bo > 8B, calculate 
4 2 
| 5/2) | 
M,i= FB ” Bo > B. (62) 
p+ —— 
B 
If Bo in (61) is < 8, calculate 
eee 3 4. 5 P (Go/2)e 9, (63) 


It is now possible to cross plot the formulas derived in this section ad 
nauseum. We shall show only two sets of graphs. Figures 10 and 11 show 


four curves of MSE/o? in dB vs V28 where E = 2aP/N» and P is the 
transmitted power divided by the parameter No/2a. In each case we 
plot the optimized results and the suboptimized results. The optimized 
decision feedback equalizer results were evaluated from equations 
(58), (59), and (60) and from equations (61), (62), and (63) for the 
linear equalizer. The nonoptimized results are given in (48) and (53) 
respectively. In all cases H = p/8 in dB. Marked on the curves is the 
value of V¥28) where the transmitting filters cut off. We show two 
cases, 10 login# = 60 and 10 login = 100. It appears that inband 
optimization does not provide a great deal of performance enhance- 
ment. As expected, inband optimization yields more improvement in 
the linear equalization scheme than in decision feedback. 

We have also evaluated the optimized results as a function of the 
actual signal-to-noise ratio p, where the noise is measured in whatever 


— 10 LOG MSE 
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band happens to be optimum. We found insignificant differences 
between these and the suboptimized results shown in Figs. 5 through 
10. 

In concluding this section we wish to stress that in practice error 
propagation problems may negate the indicated theoretical results for 
this channel. When the MSE is large, errors will result. In addition, the 
tap gains of the feedback filter may become quite large causing those 
errors which do result to propagate. 
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APPENDIX A 
Solution of the Wiener-Hopf Equations 
We wish to solve the set of linear equations 


0 
Re= Y Mi-Sn, R= Oj Ly a ey (64) 


n=—-% 


where {Ri} and Mi-: _ Rix + Non—z (6n—k = 1, n= kis On—k _ 0, 
n #~ k) are given. 

Since {M,}2. is a correlation sequence with positive Fourier 
coefficients it is well known that it can be represented as the discrete 


convolution of a sequence {M,}°,, and a sequence {M7}>, namely 
M, = p> M}Mz_,  foralln. (65) 
Let the sequence {X,}2. be determined from 
R, = y M+X.;, alk. (66) 
Substituting (65) and (66) into (64) gives 


kaw 0 
x M;* Xi-j — x SoMa = 0. (67) 
jJ= 


n=—-o0 
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Clearly a solution of 


0 
Xe= DL SiMp4y, k& 0, (68) 


n=—e 


is also a solution of (67). 
Define the two-sided discrete Fourier transform of a sequence { Xn}2. 
by 


X00) = > X,ein?. 


Take the transform of both sides of (66) to obtain 
R(6) = M+(@)X (6). (69) 


0 
The one-sided transform ( > ) of (68) is 





X~(6) = S-(6)M- (6) (70) 
and X~(6) is obtained from (69) as 
x-@ =| anos | (71) 


where [- |_ stands for “projection to negative integers only.” To 
obtain the projection, expand [- ] in a two-sided Fourier series and 
retain only the part of the series containing negative/positive coefh- 
cients (including zero). 

Thus the desired solution is 

7 | R(8) 
0 = sem Lae | oe 

To proceed further, observe that since (6) = M+(6)M-(@) and 

M (6) = R(6) + No it is possible to calculate explicitly 








R(6) _ No 
= M v0 
law | ro+S (73) 
where yo is the de coefficient of M+(6). 

The final solution for S~(@) is therefore 


No 


TO =) FOr 


(74) 

As pointed out in the text, there are various methods available for 
calculating M+(@) from a known function M(@). We briefly outline 
one such approach. Since M(6) > 0, 0 S$ 6S 27, In M(6) may be 
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expanded in a two-sided Fourier series 
0 io) 
Ins WA(G) So: Gentes ye are, (75) 
n=—0 n=0 
Knowing the sequence {y#}2. we can get immediately 
M+(0) = exp | 3 ate**| (76) 
n=0 
and 
0 . 
M-@ = exp |X ave. 
Notice that the de term of 1/+(@) equals the de term of M~(6). 


APPENDIX B 
Evaluation of Integrals 
B.1 Decision Feedback 


The detailed evaluation of 
3 
I=2 - In [1 + pe-¥880Jdy (77) 
0 


is accomplished as follows: Change the variable of integration to 
x = (vV46y — In p) which gives 


1 V28—In p 
= B Inf1 + e-*][a + In p Jd. (78) 
—In p 


Assume V¥26 > In p > 0 and write (78) as 
1 0 V2B—1n p 
: ° B Che al ) 
= = Pn p = x)de +5 f°" dnp = 2) In [1 + e7? Jdx 
0 B Jo 
1 y28—In p 
ae I (x + Inp)In[1 + e-*]dx (79) 
B Jo 
— (net (no 1 pir ge 23 
38 38 zal (In p — x) In [1 + e-* |dx 


1 26—In p 
2 I (x + In p) In (1 + e7*)da. 
B Jo 
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Since x > 0 in the range of integration, expand 
Inf tewJ= D (— Ho 
n=1 


and substitute into (79) to obtain 





as oe + re ey (— 1)" {SPA (In p) — a p) 


4, Bu(vB6 — Inp) , Inp An(V28 ~ In p) 


- ; ~ , (80) 


where 
Aa(t) = = (1 — e*8) 


and 
1 — e-™® — née~é 
n . 


Collecting terms and recognizing that 


= 1 1 
— ntl = __ 





we finally get 
_ (np)? + vIn p p)? a wr In p “vid x (— ])nH (pe Bil 


1 
n = ~ (pe-v28)n]. 
+5 ¥ (- wit (ve-8)" |. (81) 
When In p > ¥26 > 0, ae can be expressed in the form 
rs =f ‘ ee Tae) el + e* |dx 
1 s—(in e—v28) 
+ = I (a + In p) In [1 + e77 Jdx 
B Jo 


= =5f h(a 9 Sid pe ein ae geey ae 
1 In p—vV28 
ae (In p — x)[w + In (1 + e7*) Jdx. (82) 


At this stage In (1 + e-*) is again expanded in a power series and when 
the terms are collected we obtain 


T= inp — 298 + 42 ¥ (- 1)rH 


y (- 5 L — e083) 


n=l1 





a 


+ 


Wir 
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B.2 Linear Equation 


The integral we wish to evaluate here is 
3 
[T=2 I (1 + pe-Vilv)—1dy, (84) 
0 


We follow the identical procedure as in the previous case. First change 
the variable of integration to obtain 





V26—In 
r= . ” (FES) ae. (85) 
ee 


Assume that V2 > In p > 0 and write 


c A, ef: ) 


V28—In © (g + In p) Ine e~2(In p — 2) 
a Teed +5 f “(res 


_ . (v28 — In p) (28 + In p) 
+35 (— 1)*[B,(¥26 — In p) + B,(In p)] 


_— 
~ 


np 


ae (— 1)*[A,(V26 — In p) — Az(In p)] 


Ms 


+ 


Il 
ar 





=1-— es — 3+ eae 


+5 EDS (St er). (86) 
n=1 n p 
When In p > V¥28 > 0 we get 


1 si’ e-z(In p — 2) 1 in e—v28 e-2(In p — zx) 
I i Se ee Rercias AEE GSS cc 
er eer at ito * 


RBi— BDI 
Ms 


(— 1)"*"{In pA, (In p) — B,(In p) 
+ In pA,(In p — ¥28) + B,(In p — V28)} (87) 


and after collecting terms we finally obtain 


D ev28 le i i | ev28 ‘ 
[T= iat a pS be RD) ep 3 
gm [i+ p 1435! p ala (= ) 


1 


n 
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Traffic Measurement Biases Induced by 
Partial Sampling 


By A. DESCLOUX 
(Manuscript received March 22, 1973) 


Under equilibrium conditions, the sample average of the delays en- 
countered by all the calls submitted during a given time interval is an 
unbiased estimate of the mean of the delay distribution. If some of the 
delays are not observed, the resulting sample average need no longer be 
an unbiased estimator of the corresponding population mean. This is the 
case when, for instance, only a limited number of delays can be timed 
simultaneously. The purpose of this paper is to investigate these biases 
for queuing systems when only one clock is available and thus one delay 
only can be measured at a time. It is shown that, regardless of the order 
of service, the expected value of the observed average delays is always 
smaller than the mean waiting time for all calls. 

Although the average delay on all calls is independent of the order of 
service, the measurement biases resulting when only one delay can be 
measured at once depend on the queue discipline. In particular, we shall 
show that the average delay for all calls is always larger than the average 
delay of the observed calls even if these calls are always served last (ob- 
served-call served-last). 


I. INTRODUCTION 


The following remarks due to J. F. C. Kingman appear in the 
Proceedings of the Symposium on Congestion Theory held at the Uni- 
versity of North Carolina in 1964 (Ref. 1, pp. 8314-815) : ‘‘To illustrate 
the pitfalls of inference from congestion systems, let me tell a (more 
or less true) story. It was desired to estimate the mean waiting time 
in a particular queuing system, and for technical reasons, only one 
customer could be timed at once. Thus the waiting time w; of a customer 
was measured. When he entered service, the next customer to arrive 
was observed and his waiting time we was noted. This procedure 
continued, the waiting times w3, ws, --: being measured, and, for 
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large n, 


NM (wy + we + -++ + wn) 


was used as an estimate of the waiting time. It is, however, strongly 
biased and inconsistent, because of the selection of the customers to be 
observed. The mean waiting time is overestimated by a factor which 
becomes arbitrarily large as the traffic intensity approaches one.’”’ We 
stress that, according to this sampling procedure, a customer is observed 
if and only if it arrives when the clock is free. 

Another instance of biases induced by the measurement procedure 
is reported by Oberer and Riesz.? These authors have investigated the 
possibility of estimating blocking probabilities in telephone networks 
by means of test calls generated by a single source repeatedly calling 
a dedicated number. Their study shows that the proportion of blocked 
test-calls does not yield a suitable estimate of the grade of service as 
it is markedly biased downwards. As expected the bias becomes larger 
as the intervals between consecutive (nonoverlapping) test-calls 
becomes smaller. It is also established in Ref. 2 that the relative 
test-call biases increase as the blocking probability decreases. 

It is worth noting that the biases studied in Ref. 2, as well as here, 
are of a different sign than those referred to by Kingman. Much more 
important, however, is the fact that measurement techniques which, 
superficially, appear to be adequate may prove to be very unreliable. 
It is thus becoming increasingly clear that great care is required in the 
design of performance measurements for stochastic service systems 
so that unanticipated biases are not encountered. 

The purpose of this paper is to investigate the effects of partial 
sampling on the estimate of the mean (overall) waiting time obtained 
by averaging measured delays. The biases induced by such limitations 
will be studied here for M/G/1 and GI/M/s when at most one call 
can be observed at once and the estimation procedure is as described 
by Kingman. We shall see that, in these systems, the equilibrium 
average delay of the observed calls is always smaller than the equi- 
librium average delay for all calls. 

It is well known that the average delay for all calls is the same for 
all queue disciplines which are independent of the lengths of the 
individual calls (no other type of queue disciplines will be considered 
here). As we shall see, this is not true of the mean measured delay when 
only one delay can be recorded at a time. In this case, both the un- 
conditional and the conditional average delayst are (as expected) 


t As customary, unconditional and conditional delays pertain to arbitrary and 
delayed calls respectively. 
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smallest when the observed calls are served first and largest when they 
are served last. (In particular, the second extreme case occurs in 
systems with first-come last-served queue discipline.) Furthermore, 
in view of the general inequality mentioned at the end of the preceding 
paragraph, the upper bound for the unconditional average delay of 
the observed calls (which is reached when the observed calls are served 
last) is always a strict lower bound for the unconditional average delay 
for all calls! 

The preceding result pertains to unconditional delays and does not 
always hold for the average delay of those sampled calls which en- 
counter a delay. Thus for M/M/s, the conditional average delay for 
all delayed calls is equal to the conditional average delay of the ob- 
served delayed calls so long as these are always served last. (Note that 
for M/M/s, the average delay of the delayed calls is equal to the 
average length of the busy period, and that the waiting-time distri- 
bution of the observed calls coincides with the busy-period distribution 
for the observed-served-last measurement procedure.) In contrast, 
for the M/T;,/1 queue, the upper bound for the conditional average 
delay of the observed delay calls is larger than the average conditional 
delay for all delayed calls when k > 1, the inequality being reversed 
whenever 0 < k < 1. (1; is used here to designate the gamma distri- 
bution with mean 1 and variance k~!. Thus M/T;/s is identical to 
M/M/s. When &k is an integer, I, is the Erlangian distribution often 
designated Ey.) 

Expressions for the moments of the equilibrium delay-distribution of 
the observed calls are given for M/G/1 and first-come first-served 
queue discipline. The equilibrium delay-distribution of the observed 
calls is also derived for M/M/s with order-of-arrival service. (Cor- 
responding results for the ‘‘observed-call served-first’’ and the ‘‘ob- 
served-call served-last’’ measurement procedures are immediate.) These 
formulas are used to show that the biases induced by partial sampling 
can be quite substantial. 

When the average service-time is unity, an assumption made 
throughout, the average delay, EW, for the single-server queue 
M/G/1 is given by the formula (Ref. 3, pp. 46-50) : 


EW = am2/2(1 — a), 


where a is the server occupancy and mz is the second moment about 0 
of the service-time distribution. Since mz can be arbitrarily large, no 
bound can be placed on the value of HW. But when only one delay 
can be timed at once, we shall see that the expectation of the observed 
delays cannot exceed 1/(1 — a). Therefore for any prescribed value 


1378 THE BELL SYSTEM TECHNICAL JOURNAL, OCTOBER 1973 


of a, it is always possible to find service-time distributions for which 
the ratio of the average delay for all calls to the average delay of the 
observed calls exceeds any given bound. 

To simplify the exposition we restrict ourselves to full-access delay 
systems with recurrent inputs in which delays are measured by means 
of a single clock. Some of the results obtained below can, however, 
be extended to more general situations. 

(In the sequel, W, with or without affix, is used to designate the 
waiting time of an arbitrary call while W,, with or without affix, is 
used as the generic symbol for the observed delays when only one 
clock is available.) 


II. A GENERAL DELAY FORMULA 


Consider the queuing system GI/G/s and suppose that the arrival 
and service-time distributions are such that equilibrium can be 
reached. (To avoid trivial qualifications we assume throughout that 
the mean interarrival time is finite. For the same reason, the under- 
lying distributions are also supposed to be such that simultaneous 
occurrences of events need not be considered.) The purpose of this 
section is to derive a formula relating the average delay of the observed 
calls to the equilibrium probability, 6, that an observed call has 
immediate access to a server. To this end we prove first that 


eSB) A)/ LO SB) A) Bl, (1) 


where B is the equilibrium blocking probability for all calls and A is 
the expectation of the number of unobserved calls originating during 
the waiting time of an arbitrary observed delayed call. 

It follows from (1) that the probability that an observed call is 
blocked is always (strictly) smaller than the overall probability of 
delay (so long as B ~ 0 or B ¥ 1, two trivial cases that we exclude 
from our considerations). This, of course, is a consequence of the fact 
that all nonblocked calls are observed whereas, with one clock only, 
some delays may not be recorded. 

We turn now to the proof of (1). Consider an infinite sequence of 
consecutive calls and for the 7th call (¢ = 1, 2, ---) let 


ie 0 if the zth call is delayed, 
: 1 if the zth call is not delayed, 
0 if the zth call is not observed, 
Y; = <0 if the zth call is observed and not delayed, 


1 if the 7th call is observed and delayed. 
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Let « > 0. Then assuming that the system is in equilibrium when 
the first call arrives, we have, by the integral stationarity theorem 
(Ref. 4, p. 419), 





iy Ma Spee Rok EX, (2) 
noo (Xi + Yi +e) +-- + (X, + Ys Fe E(Xit Yi t+. 


and 








n70 (Xi+te +---+ (X,+ ¢ HX, + « , 


with probability 1. However, since (Ref. 4, p. 421) 


lim “2+ 7 0. px, = Pr[X, = 1] > 0, 


with probability 1, there is, for almost all realizations of the process, 
an integer n such that the ratios 


Xit-::+ Xn ene: 
(Xi + Ys) +--+ (Xm + Ym)’ wee 


are well defined. Hence, by (2) and (3), we have 


EX, lim Xit-:::+ Xa 
1 See oi ee eee (Xi + Yi) +:: -+ (X, + Y,) 
< EX, + € 
= H(Xi+ Vi)’ 
with probability 1 and, letting « tend to 0, 


i oe aes EX, 
nove (Xa Ya) +--+ (Xn + Yn) EX FY)’ 


with probability 1. 

[The preceding derivation makes use of the fact that the station- 
arity—and hence the integral stationarity—of the processes {X,, 
a= 1,---} and {X;+ Yi,7 = 1, ---} follows from the property that 
the random variables X; and Y;,7 = 1, ---, whose means are finite, are 
“translates”? defined on the stationarity queuing process (Ref. 4, 
p. 417 ff.). 

[The formulas in Ref. 4, p. 419, Theorem A, involve conditional 
expectations with respect to fields of invariant events. Under the 
present circumstances these expressions can be simplified. Indeed 
let us specify the state of the system, 7; at time t, by means of the 
vector whose components are the arrival time of the last request 
placed before ¢t and the elapsed portions of the service-times in progress 





(4) 
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at time ¢. Then the only invariant sets (Refs. 4 and 5) of the process 
{J;, — © <t< ©} are the whole space and the null-set. This 
property, in turn, implies that the conditional expectations of the 
random variables X; and X;+ Y; relative to their invariant fields 
can be replaced by the unconditional expectations HX, and H(X; + Y,), 
respectively. These and other similar substitutions are made here 
without formal justification. ] 
Consider now an infinite sequence of observed calls and let 


Gites 0 if the 7th observed call is delayed, 
: 1 if the 7th observed call is not delayed. 


Then (Ref. 4, p. 421) 


Zit Zot---+Zn 
n 


nan 


Since (4) and (5) are both equal to the proportion of observed calls 
with zero delay over an interval of infinite length, we have 


We note that EX, is equal to the probability that a call (observed 
or not) is not delayed. Hence 


EX,=1-B, (7) 
and to complete the proof of (1) we have to show that 
EY, = B/(1+ A). (8) 


To this end consider again a stationary sequence of observed calls and 
let A; be the number of unobserved calls placed during the waiting 
time of the 7th call that is both observed and delayed. Then we have 
(Ref. 4, pp. 419-421) 


Lit ae A Pec Ag Tea (9) 


with probability 1. 
Furthermore, by the integral stationarity theorem (Ref. 4, p. 419) 
and a simple e-argument of the type used in the proof of (2), we have: 


faa Va ise Yin EY, 
pos: ee Xi) tee + (1 X n) ~ EG — Xi)’ 


with probability 1. 


(10) 
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Since the left-hand sides of (9) and (10) are both equal to the propor- 
tion of delayed calls that are observed, we have 


EY, = E(1 — X})/(1+ A) = B/(1+ A). 


This completes the proof of (1). 

Our next step will now be to relate A to the average delay, EW,,, 
of the observed calls. Let W;, be the delay of the 7th observed call 
and let U; be the interval between the end of the 7th and the beginning 
of the (¢ + 1)st measurement. Let also J, be the interval between 
the arrival epochs of the nth and (n + 1)st call (in the whole sequence 
of calls, observed or not). Then we have: 


(Wig + Ui) +++ (Way + Un) = Ti +e+ ++ Ie, (1) 


where Kz, a random variable, is equal to the number of calls placed 
during the interval that starts with an observed call and ends just 
before the beginning of the (n + 1)st measurement. By the stationarity 
theorem, we have (Ref. 4, p. 421): 


Him (Wie + Uy) Fe (Wg + Un) _ E(Wi,+ Uy), (12) 


pee n 


with probability 1, and, by the strong law of large numbers (note 
that the J,’s are, by assumption, independent random variables with 
finite means and that K, = n), 


ive dpa cesar Te 
n> eo K, 


with probability 1, where a—! is the expected interarrival-time. 
Furthermore, 


so that 


lim “* =EZi+0-Z)1+4A)] =8+ 1-811 +4), 14 


=o, (13) 


with probability 1. 
Combining (11)—(14) we find that 


aE (Wiz + Us) =1+A(1 — 8). (15) 


In particular, when the input is Poissonian, HU; = a and (15) 


reduces to 
alWi, = A(1 — ®). (16) 
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Thus, taking (1) into account, we find that: 
EW, = EW,, = (®+ B— 1)/a(1 — B). (17) 


It should be noted that the preceding relation is valid regardless of 
the order of service. 

When the calls are served in order of arrival, (16) is an immediate 
consequence of the fact that the waiting time of any given call is not 
affected by the stream of requests placed after its arrival epoch. This 
is also true when the observed calls are always served first and (16) 
can then be written down with equal ease. 

We note that (17) can be obtained quickly whenever the epochs 
at which measurements begin or terminate constitute a renewal 
process. In such cases, the expected number of observations in time ¢ 
is (asymptotically) 


(EW, + a)-!-¢ + 0(1), ¢ large, (18) 


and the expected number of arrival points at which the system is 
empty in time ¢ is 


a(1 — B)-¢ + 0(1), ¢ large. (19) 


The long-term proportion of observed calls with no delay is given by 
the ratio of (19) to (18) with probability 1 (cf. Ref. 6, p. 264, alter- 
native form of Theorem IV): 


a(1 — B)(EW, +o). (20) 


Since the probability @ that an arbitrary call is not delayed is in- 
dependent of the past, a simple application of the strong law of large 
numbers show that (20) may be equated to ® and (17) therefore holds. 
An instance where the preceding argument can be applied is the 
M/G/1 system with observed calls always served last. With this 
order of service, there is exactly one call in the system at the termina- 
tion of each measurement. These epochs constitute a renewal process 
since they also coincide with the beginnings of the service-times of 
the observed calls. With an obvious change, the previous argument 
remains true for M/M/s with observed-calls served-last. 


II. TWO EXTREME CASES 


In this section we show that if only one clock is available then the 
expected average delay of the observed calls is largest when the 
observed calls are served last and smallest when the observed calls 
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are served first. These relations are not statistical: They are satisfied 
by all the realizations of the process over any finite or infinite time 
interval regardless of the arrival and service-time distributions. (To 
avoid ambiguities, we assume that the timing device is free at the 
beginning of the realizations.) 

We note first that under any measurement procedure, all the calls 
which arrive when all the servers are busy, but no request is waiting, 
are observed. These calls are the only delayed calls that are observed 
when the observed calls are served last. Therefore (7) the number of 
observed delayed calls takes its smallest value for the observed-served- 
last procedure and (27) during the measurement of a delay, D, under 
this particular procedure any observed delay under any alternate 
single-clock measurement procedure cannot exceed D. Combining these 
two facts and taking into account that all calls with zero delay are 
observed we may conclude that the observed average delay takes 
always its largest value when the observed calls are served last, as is 
the case for the first-come last-served queue discipline. 

When the observed calls are served first, we note that (7) the number 
of observed delays over any busy period (initial or not) is never smaller 
than for any other single-clock measurement procedure and (iz) to 
each observed delayed call there corresponds, under any other single- 
clock measurement procedure, one call whose delay is at least as 
large and this correspondence involves all the observed delays under 
the alternate procedure. All calls with zero delay are again observed 
and the average delay of the observed calls takes therefore its smallest 
value when the observed calls are served first. 

Clearly the conditional average delays of the observed calls do have 
the same property. 


Iv. BOUNDS FOR THE AVERAGE DELAY OF THE OBSERVED CALLS IN 
M/G/1 


The object of this section is to determine the upper and lower 
bounds for the average delay, EW,,, of the observed calls in M/G/1. 
These bounds, as we have seen, are reached when the observed calls 
are served last and first respectively. Under the present conditions, 
formula (17) may be written as follows: 


EW, = (®+a—1)/a(1 — a). (21) 


For a given server occupancy a, HW, is a monotone increasing function 
of @ = ®(a). Since ® < 1, (21) implies that EW, < (1 — a). Hence 
fora < 1, EW,, is always bounded (but EW is not). 
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It will be convenient to define here the service backlog, at a given 
instant ¢, as the sum of the service-times of all waiting requests plus 
the residual of the service-time of the request being served. [When 
calls are served in order of arrival, the service backlog is equal to the 
virtual waiting time (Ref. 3, p. 59 ff.).] 

Now let F(-) be the stationary cumulative distribution of the 
service backlog at the end of a measurement. The probability, (a), 
that an observed call does not suffer a delay is simply the Laplace- 
Stieltjes transform of F(-) evaluated at a since it is equal to the proba- 
bility that no call originates during a time interval whose length is 
that of the service backlog: 


b(a) = i ” etd F(t), 


Writing o(-) for the Laplace-Stieltjes transform of the service-time 
distribution we have the following inequality : 


P(a) S o(a). (22) 


This inequality is a consequence of the fact that, at the conclusion of a 
measurement, the service backlog may be represented as the sum of 
two independent random variables, one of which is the full service-time 
of the request whose delay has just come to an end while the other is 
equal to the sum of the service-times of all the waiting requests. 
Writing R(-) for the c.d.f. of the latter and S(-) for the service-time 
distribution, we have: 


(a) 


I 


T etd [RU ~ dS) 


o(a) i ” e-*tdR(t) S o(a). 
0 
When the observed calls are served last, 


1 for t2 0, 
a 10 for t <0, 


and 
@(a) = a(a). 


We can therefore conclude that 


ote) eet, 


EW, 3 (23) 
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We are now in a position to prove that the average delay, EW,, is 
always smaller than the average delay for all calls (observed or not). 
Since the service-time is unity, we have: 





a?M2 
’ 


bi) = Gays os emdS() <1-at 3 


where mm: is the second moment of S about the origin. Hence, substi- 
tuting 1 — a + a?m,/2 for o(a) in (23) we find that, irrespective of the 
service order: 


ame ee 
EW, < xy = EW. (24) 


By (23) and the equality in (24) we also have: 


EW/EW, > see 
so that, for any given a, we can always find a service-time distribution 
such that ELW/EW,, exceeds any preassigned value. 

We now derive an absolute lower bound for EW,,. As shown above, 
this bound can be found by assuming that the observed calls are 
served first. Our first step here is to determine A. 

For the observed-served-first procedure, the circumstances under 
which a positive delay can be observed are as follows: at some time 
the clock is not in use and a service-time begins, and during this service- 
time a new call arrives. Thus A is the conditional expectation of the 
number of arrivals minus 1 during an arbitrary service-time given that 
at least one call is placed during a service-time. A, therefore, is given 
by the formula 


ioe] ie} 


A 


(n — 122" --otas (i) / [ (1 — e-**)dS(t) 


0 n= 


i ” Cot — 1 + e-**Jas(2) Vi i * 1 — e*)dS(t) 


[oe — 1+ o(a)]/L1 — o(@)]. 


By means of (1) and (21), it is now readily shown that, for the ob- 
served-calls-served-first procedure : 


ll 
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and 
_ lay a 1 
se ae eormanerr eye 
Summing up, we have the following inequalities for @ and EW,, 
regardless of the measurement procedure: 


l—e 
tea ey a a 
and 
_o(a) Fa = 1 ee) ie 
a[2 — a — o(a)] Bet a(1 — a) ) 


For exponential service-times, (25) and (26) reduce to 


(1 —o’)/1 +a—o’) SO (a) S1/(1 +a), 
a/(1+a—o’) S EWys S a/(1 — a’). 


(The subscript 1 is added to ® and EW,, to indicate that the service- 
times are exponentially distributed.) 


V. THE SINGLE-SERVER QUEUE M/G/1 WITH ORDER-OF—ARRIVAL 
SERVICE 


In this section we consider the M/G/1 queue under the assumption 
that the calls are served in order of arrival. Our principal aim here is 
to determine 6 = ®(a) and then, by means of (17), EW,,. To this end 
let pn be the probability that there are n calls in the system immediately 
after the conclusion of a measurement. Then, relating the state proba- 
bilities at two consecutive conclusions of delay measurements, we 
find that 


Di = > p» fi e—@td § (™) (¢) + x po fe ate~*"dS ™ (£), 
n=1 0 n=l 0 
0 ee) k 27 
Dr = PS ps | cy etd S(™) (£), k> 1, ( ) 
n=1 0 . 


where S“ is the nth convolution of the service-time distribution, S, 
with itself. 
Now let 


G(z2) = D pa 


Equations (27) yield: 
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G (2) Dm&x™ = 3 “™ 3 p» | (at)"_ 2148) (t) 
m=1 n=1 o mM. 


a5 We i * e-atd im) (2) 
n=l1 0 


il 
iMs 


1 


I 


af. | oem asm + 2 E pao) 


m=1 


Me iMs 


Dn / : eer? = Dds Gy) she S pro” (a) 
0 n=1 


ll 
— 


n 


& pootfa(l — 2)] — (1-2) ¥ pro*(a) 
G{oLa(1 — x) ]} — (1 — 2)G[o(@)]. 
Summing up, we have the relation 
G(x) = G{o[a(1 — z)]} — (1 — 2)G[o(@)]. (28) 
Note also that 
Bla) = psf” eS) = Glo(a)] 


Let 20 = o(a) and z, = o[a(1 — 2,1) ],n = 1, 2,---. 
Since 0 Sa <1, we have rp < 41 <---S 1 and lim,..2, does 
therefore exist. With this notation, we obtain, from (28): 


B(a) = G(a1) — (1 — 20) PQ), 


G(x) = G(x2) — (1 — 21)®(@), (29) 
(ays) = Glts) —  — 24-4) 8(0). 
Adding up these relations, we find that: 
#@)|1 +E a — 2m) | = C(es), 
and, by passing to the limit, 
(a) | 1 is ie em) | 4, (30) 


(Note that lim,..G(x,) exists since the x, constitute a positive 
monotone-increasing sequence bounded by 1. By letting n —© in the 
last of the relations (29) it follows immediately that limnz..%2 = 1 so 
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long as ®(a) #0. This last condition is however clearly satisfied 
whenever a < 1.) 

In particular, when the service-times are negative exponential, we 
have: S(é) = 1 — e', t = 0, c(s) = (1 +s)", and (1 — tm) = a™t!/ 
(1 +a+---+ a+), Hence, by (80), 


® eo een a 31 
@) = | ee ee el) 

We examine briefly the case where the service-times have a gamma 
distribution with parameter k (the subscript k is added to the symbols 
considered earlier in order to stress their dependence on k). We have, 
in this case: 


ox(a) = [k/(k + @) I. 


Then oz(a@) is a strictly decreasing function of k(> 0) as can be seen 
by taking the derivative of In of !(a) = In (1 + a/k)* and using the 
inequality In (1 — 2) >2/(1+2), «>0 (Ref. 7, p. 68). This 
monotonicity property of o;, implies that 


LKo > Leth, Lia Leth ay 5 k>0, h>0, 


and we have, therefore: 


>, (1 as Lin) < 2d, (1 = Lkth,m) h> 0, 


m=0 
so that 
Bi(a) > Piin(a), 
and from (21) 
EW yx > EW y,ntns k > 0. 


For k = 1 (exponential service-time) the conditional average delay 
for the delayed calls under the observed-last-served procedure is equal 
to the average length of the busy period. To see this one need only 
note that each positive observed delay begins with an arrival that 
occurs when there is exactly one customer in the system and ends 
when, for the first time thereafter, there is no waiting customer. 
Hence, for k = 1, the conditional delay distribution of the observed 
calls is the same as the busy-period distribution (Ref. 8, p. 32). Since 
the average length of the busy period and the average of the condi- 
tional waiting times of all the delayed calls are both equal to (1 — @)~! 
(Ref. 3, p. 63), we have 


PU BMS NO ey, GR) 


TRAFFIC MEASUREMENT BIASES 1389 


where HW,.1 designates the average delay when the observed calls 
are served last. 
Clearly, the inequalities (22) and (23) imply that 


LEW xx o.(a) ta — 1 
1 — ®, ~ [1 — o;(a@) Ja(1 — a) 


We note that if the service-times have a gamma distribution with 
transform o;(s) = (k/k + s)*, then the conditional average delay on 
all delayed calls is given by the formula (Ref. 3, p. 50): 

EW, _[k+1 1 
a k 2(1 — a) 
Substituting (k/k + a@)* for ox(a) in (83) we obtain the following 
upper bound for the conditional average delay of the observed delayed 
calls (this bound, as we know, is reached when the observed calls are 

served last) : 








(33) 





(34) 


[k/kta)eta-1l 
{1 — [k/(k + a) }¥}a(1 — a) 
Subtracting (34) from (35) we find that the difference is of the same 
sign as 


(35) 





a— {1 —[k/(k +a) }}(1 + (+ Io/2]. (36) 


The two factors in the second term of (36) are both increasing func- 
tions of k and since (36) vanishes for k = 1 [a fact that we already 
know from (382) ] we may conclude that (35) is smaller than (34) for 
0<k <1, and greater than (34) for k > 1. This proves that, for 
k < 1, the conditional average delay for all delayed calls is still larger 
than the conditional average delay for all observed delayed calls even if 
these calls are served last. For k > 1, the conditional average delay 
of the observed delayed calls for the observed-served-last procedure 
is larger than the conditional average delay of all the delayed calls. 

Expressions for the higher moments of the observed calls delay- 
distribution are also readily obtained. Let F be the equilibrium cumu- 
lative distribution of the virtual waiting time, V, at the conclusion 
of the measurement of a delay (this delay, of course, may be equal to 
zero). The delay distribution, K, of the calls whose delays are observed, 
can be readily expressed in terms of F’. Indeed we have: 


K(w) = PrOW, sw} =af”| [exp — aly ~1)-dr lat 


a: i ” e-audF(y). (37) 
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TaBLE I—MEANS AND STANDARD DEVIATIONS OF THE DeLay Dis- 
TRIBUTIONS FOR ALL CALLS AND FOR ALL OBSERVED 
Cauts (1 Ciock) In M/M/1—Frrst- 
ComE First-SERVED 


[e4 EW, SWi EW, SW 
0.1 0.11111 0.48432 0.09259 0.42264 
0.2 0.25000 0.75000 0.17840 0.58463 
0.3 0.42857 1.0202 0.26684 0.72360 
0.4 0.66667 1.3333 0.36669 0.87308 
0.5 1.0000 1.7321 0.48958 1.0590 
0.6 1.5000 2.2913 0.65554 1.3188 
0.7 2.3333 3.1798 0.90754 1.7310 
0.8 4.0000 4.8990 1.3659 2.5191 
0.9 9.0000 9.9499 2.5808 4.7519 
0.92 11.500 12.460 3.1398 5.8266 
0.94 15.667 16.637 4.0284 7.5793 
0.96 24.000 24.980 5.6974 10.987 
0.98 49.000 49.990 10.243 20.776 
0.99 99.000 99.995 18.393 39.434 


The problem of finding the distribution K of the observed delays 
is thus reduced to the problem of finding F. The distribution F satisfies 
the following integral equation: 


PO =% fam f- gay OD” apy) 
ay i: * dS (u) f ° ev) (y), (38) 


where S“ stands for the nth convolution of S with itself. Equation 
(88) follows immediately upon noticing that, at the conclusion of a 
measurement, the only calls in the system are (2) the call whose delay 
has just been measured and (77) all the calls which arrived during the 
measurement interval (we note that the preceding argument makes 
essential use of the assumption that calls are served in order of arrival). 

Let ¢ and o be respectively the Laplace-Stieltjes transform of F 
and S. Then transforming (38) we obtain 


(8) = ¥ or(s) f° w P ariy) + o(@) f° mary) 
= {all — o(s)]} + S@)[o(s) — 1]. (39) 
This formula may be used to derive recurrence relations for the 
moments of V. Let niu, = EV” and 
d 


Nn = (— 1)" a 7 (8) 





= I ” tdS(t), 
0 0 


s= 
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Fig. 1—Average delay for M/M/1 vs occupancy. 


so that n!v, is the nth moment of the service-time distribution. Using 
Faa di Bruno’s formula for the derivative of a composite function 


(Ref. 9, p. 36) we find that: 


k! 
bn = s kyl --k, Tee Lee ‘ a =F D(a) Yn, 


(40) 


with k = ki +---+ k, and the sum over all solutions in non-negative 
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Fig. 2—Conditional average delay for M/M/1 vs occupancy. 


integers of ki + 2k, +---+ nk, =n (note that »v; = 1 since the 
average service-time is assumed here to be equal to 1). 
In particular, for n = 1, 2, and 3 we have: 


11 = BV = 8(@)/(1 — a), 
a = Sr = ®)ne/(1 — @)(1 — 


ys = = = @(a)[BaXv} + va(1 — a) ]/ (1 — a) (1 — a) (A — a) 


where (qa) is the probability that an observed call is not delayed. 
When the service-time distribution is exponential (with mean 1) 
we have »; = 1,7 = 0,1, ---, and (40) becomes: 
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Fig. 3—Average delay for M/D/1 vs occupancy. 


k! 
aoe De iis kt Bee + #;(a) 


pa ae 1) meat + (2), pea: 


Equation (87) may be used to express the moments of K in terms of 
the moments of V. We have: 


EW? = i w'dK(w) = af” wr [" exp — a(y — w)dF (y)du, 
0 0 w 
and, upon integrating by parts, we obtain 


EW, = EV - ae ~ B(a)], 
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Fig. 4—Conditional average delay for M/D/1 vs occupancy. 





and 
Ewe = Even "tly ado. 
Qa 
Thus, in particular, we have 
_ €@) +a-1 
UL a a(l—a) ’ 


a Cea) Saltese = 1 Se) 
EW. a(1 — a)(1 — a) 


The moments of W,, depend only on the moments of the service-time 
distribution and on ®(a). 

As a numerical illustration of the biases induced when only one 
clock is available, the means and the standard deviations of Wi and 
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Fig. 5—Conditional average delay for M/T;/1 vs occupancy. 


W,.1 are given in Table I. (The standard deviations of Wi and Wy 
are designated by SW, and SW,, respectively.) For further quanti- 
tative results, see Figs. 1-6. 


VI. THE SINGLE-SERVER QUEUE M/M/1 


In this section we consider a single-server delay system and assume 
that: (¢) calls arrive in a Poisson process of intensity a; (a) the 
service-times are independent random variables with the same negative 
exponential distribution; and (diz) calls are served in order of arrival. 
We again suppose that only one delay can be measured at a time. Our 
purpose here is to derive the delay distribution of the observed calls. 

Let F(-) be the equilibrium cumulative distribution of the virtual 
waiting time at the conclusion of the measurement of a delay; the 
measured delay may of course be equal to zero. From (88), the distri- 
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Fig. 6—Conditional delay distributions for M/M/1 — a = 0.5 first-come first-served. 


bution F(-) satisfies the following integral equation: 


F(t) = a | ‘ = enw a CO du 


ns 





é [ ‘ e-udu- f ° eau F (y). 
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This relation implies that F(-) is continuous and that the virtual 
waiting time, at the conclusion of a measurement, has a density func- 
tion f(-) Lat ¢ = 0, the latter is defined as the right-hand derivative 
of F(-) ]. We have therefore 


= / : flyje™ = ca)" — ydy + et i} ; f(y)e-edy 


n! (n 


Il 








F (6) 


abet [* f(y)ererytLiL2(aut)*ldy + ef fy)edy, (41) 
where J/;(-) is the modified Bessel function of order 1 (Ref. 7, p. 374). 
The preceding relation implies that f(-) is of the form 


f® = fad) + ces, 


where c is a constant. Substitution of this expression in (41) yields, on 
taking relation 29.3.81, p. 1026, of Ref. 7 into account: 





full) = abet [ faye ALB ay)" Vdy + Ge erro. (42) 


Thus fi(-) is of the form 
filt) = Fa(t) ae ‘al + a)? exp t/(1 a a), 
and substituting this expression in (42) we find that f2(-) is of the form 


fot) = fa(t) + exp — t/(1 +a + a”). 


eae, 


Proceeding in this manner, we define successively f1(-), fs(-), --+, and 
it is readily shown, by induction, that: 


ca™ 


Fm(t) = fmti(t) + Pata oso 
‘exp — t(//(lt+a+t+-:---+ a”), m= 0,1, 2, ---; 
fo(-) = f(-). (48) 
Passing to the limit, we obtain, in this manner: 


a™ 


f®) = fo(t) +e ps fee ery 
‘exp — t//(l ta +--:+ a”), 


where fo(+) = limmsofm(-) satisfies the integral equation: 


fo(t) = abet [© fa (yereryiTs[2 (ayt)* dy. (44) 
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Note that, by virtue of (48), fm(t) > fm4i(t) for all t and that 
lim. m+eofm(t) does therefore exist and is non-negative since fm > 0 for 
all m. 

We shall prove now that the only non-negative solution of (44) 
is fo(t) = 0. 

Let 


a(s) = i ” 5, (ten*dt. 


Then, transforming the previous relation, we have by Ref. 7, p. 1026, 
equation 29.3.81: 


I 


6(s) at f° eres [ fa (yey? Ti 2 (ayt)? |dy - dt 
0 0 


= af foo (y)e-euyh (xy) 4 (eel — 1) dy 
0 


i fo(y) exp — ay (1 ae) i ay 


= I fo(y) exp (— ay)dy 


6 (=) — 6(a). 


Setting s equal to zero in the preceding relation, we find that 6(a) = 0 
which implies that f..(t) = 0 and we have, therefore, for exponential 
service-times : 





HO «portance, 
-exp — (//1 +a+---+ a”), t 
2 = a (45) 
a ea p2 oltat-::-+a™ 
‘exp — t(/(l t+a+:--+ a”), t 
f@) = F@® = 0, t<0, 
where the constant c is determined by the condition F'(0) = 0. 
By means of (37) and (45), it is readily seen that the conditional 
delay distribution is given by the following formula: 


IV 
= 


IV 
= 


qmtl 
ee ee 
-exp —t//(l1ta+--:-+a™), t>0, 


Pr [W,, = t|observed call delayed] = c’ ps 
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where c’ is determined by the requirement that the preceding expres- 
sion be equal to 1 for t = 0. 

The effect of partial sampling on the delay distribution is illustrated 
in Fig. 6. 


VII. THE MULTISERVER QUEUE M/M/s 


In this section we consider a full-access multiserver delay system 
with Poisson arrivals and exponential service-times. Our purpose here 
is to determine the probability ®{° that an observed call is served 
without delay and the expected delay HWS} of the observed calls 
(Of) = &,, WY) = W,,). This is easily done. Indeed, under the present 
assumptions, A{® the expected number of nonobserved calls during 
the measurement of a positive delay is: 


aw ya 


(s) — 
A d s(1 = ®}) , 


(46) 
where EW, and ®,; pertain to the single-server queue with load a/s. 
Equation (46) is an immediate consequence of the fact that in an 
s-server system with demand rate @ and service rate 1, the conditional 
average delay of an observed call, EHW${/(1 — ®{), is equal to the 
conditional average delay of an observed call in a single-server queue 
with offered load a/s and service rate s, i1e., HWy1/s(1 — 1). 
Hence, by (46) and (1) we have 


_ = BLG/)EWat1- ei) 
@— Bl@/)EWa+1—%]+B0 —&) 


so that, by (17), 


a 


B-EWyi 
(s) — Soe ES 
HW 5 =e) all = B)EWa 
We note that 
EW (1 — a/s)EWya 


EWP” 1—4,+ (0 — B)e/s)EWy 
For a/s fixed, the blocking probability B is strictly decreasing and 
tends to 0 as s increases. Hence 


J (s) (s-++m) 
ee eo ; m > 0, 
and 
Faia Stine oe Sees. 
EW |e sox EW 1-414 (e/s)EWyi 


1400 THE BELL SYSTEM TECHNICAL JOURNAL, OCTOBER 1973 


We stress that the preceding relations are valid regardless of the order 
of service. 

Numerical values are given in Table II. They show, in particular, 
that, for a given server occupancy, the magnitudes of the relative 
biases become larger as the number of servers, s, increases but remain 
bounded. 


VIII. AN INEQUALITY FOR GI/M/s 


For the M/G/1 queue we have seen that the average delay on all 
calls, EW, is always larger than the expected delay EW,, even if the 
observed calls are served last. It will be shown here that the same 
relation also holds for the multiserver queue GI/M/s. 

When the observed calls are served last, the waiting times of the 
observed delayed calls have the same distribution as the busy period 
whenever the service-times are exponential. Writing EW} for the 
unconditional average delay for the observed-served-last procedure 
we have therefore: 


EWS, = (1 — &)/s(1 — d), (47) 


where 6 is the root of smallest absolute value of the equation (Ref. 10, 
p. 225 ff.) 
z= A*(1 — 2) 
and A* is the Laplace-Stieltjes transform of the interarrival distri- 
bution A. We note here that 6 is also the blocking probability in the 
associated single-server queue GI/M/1 with A as interarrival distri- 
bution and exponential service-time distribution with mean 1/s. 
By means of (1) we can rewrite (47) as follows: 


EWS = B/s(1 — b)[(1 — B)(1 + A) + BY, (48) 


where B is the probability of delay in the GI/M/s queue. 

When the observed calls are served last, 1+ A is equal to the 
expected number of calls served during a busy period of GI/M/s 
which, in turn, is equal to the expected number of calls served during 
a busy period of the associated single-server queue GI/M/1. Hence, 
we have (Ref. 10, p. 286): 


1+ A= (1 — b)-, 
and on taking this relation into account, (48) reduces to 


EW = B/s[1 — B+ B(1 — b)]. 


TasLE JI—AveErRAGE Dr.uays In M/M/s—First—-ComE First-SERVED 


a/s | EW | EWu/EW,| EWS | EWR/EWS | ews | eEwSvew? | ewe | ew2yew® | CeWe/EW1). 


0.1 0.093 0.83 0.008 0.82 0.82 
0.2 0.178 0.71 0.029 0.70 0.002 0.69 0.69 
0.3 0.227 0.62 0.059 0.60 0.008 0.58 0.0004 0.57 0.57 
0.4 0.337 0.55 0.099 0.52 0.018 0.49 0.0019 0.48 0.48 
0.5 0.490 0.49 0.151 0.45 0.037 0.42 0.0059 0.40 0.39 
0.6 0.656 0.44 0.224 0.40 0.065 0.36 0.0146 0.34 0.31 
0.7 0.908 0.39 0.336 0.35 0.111 0.31 0.0316 0.28 0.24 
0.8 1.37 0.34 0.541 0.30 0.199 0.27 0.0665 0.23 0.16 
0.9 2.58 0.29 1.09 0.26 0.438 0.22 0.166 0.19 0.086 
0.95 4.71 0.25 2.06 0.22 0.866 0.19 0.349 0.17 0.045 





SHSVIG CLNAWAYOASVAW OIAVUaL 


IOFT 
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But the average delay for all calls is given by the formula (Ref. 11, 


p. 383): 
EW® = B/s(1 — b) 
so that — 
EWS, < EW?. 
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Use of a Gate to Reduce the Variance of 
Delays in Queues With Random Service 


By R. D. COLEMAN 
(Manuscript received April 18, 1973) 


We consider an N-server queuing system with Poisson arrivals and 
exponential service, in which arriving customers must pass through a 
gate into a waiting room before becoming eligible for service. Customers 
who find the gate closed wait outside until the gate opens; customers 
inside the waiting room are served at random. When the last customer 
inside acquires a server, the gate admits all those outside and then closes 
again. If no customer is waiting outside when the gate opens, the gate 
remains open until there ts a queue of k waiting customers. 

Service offered by this system is intermediary between random service 
and order-of-arrival service. As long as the gate 1s open and fewer than 
N + k customers are in the system, service is purely random. The param- 
eter k can be regarded as a threshold at which the queue is judged too long 
to permit random service to continue. 

Our main results are (1) the Laplace-Stieltjes transform of the equilib- 
rium distribution of the waiting time of an arbitrary customer and (it) 
a comparison of the second moments of the waiting time for different 
values of k with those of the waiting time under random service and order- 
of-arrival service. The service is shown to be ‘‘nearly random” at low 
loads and “not quite order-of-arrival” at high loads; for higher values of 
k this transition occurs at higher traffic intensities. 


I. INTRODUCTION 


Switching systems, particularly electromechanical switching systems, 
are often constructed so that if several customers are awaiting service 
simultaneously, they will receive service in what is essentially random 
order, i.e., a server which becomes idle will choose its next customer 
at random from the queue of customers awaiting service. Such an 
arrangement may be satisfactory when the traffic intensity is low) 
but as the intensity increases, a progressively greater number of 

1403 
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customers will have to wait an undesirably long time, and the quality 
of service may become unacceptable. 

There are available, however, methods of providing service other 
than ‘‘random service’; the most obvious one is service in order of 
arrival. The quality of service, which depends in part on the variance 
of the waiting time, will still diminish as the traffic intensity increases, 
but not as quickly as when random service is used. In fact, order-of- 
arrival service is the “best” discipline (at least when the order of 
service does not depend on individual service times) in the sense that, 
for a given traffic intensity, and hence mean waiting time, the variance 
of the waiting time is smallest.! Unfortunately, it may not be worth- 
while (or even possible) to build a system which offers service in order 
of arrival. One is therefore led to consider an intermediary queue 
discipline, one for which the variance of the equilibrium waiting time 
lies between that of random service and that of order-of-arrival 
service. 

Suppose we have an M/M/N queuing system with customers arriv- 
ing at rate \ and requiring a mean service time y~!. Suppose further 
that customers must first pass through a gate into a “corral,” or 
‘“‘waiting room,’’ before becoming eligible for service. So long as there 
are not more than N customers in the system, the gate remains open; 
an arriving customer enters the corral immediately and, if some server 
is idle, begins service. AS soon as a customer has to wait (having 
found N customers in the system), the gate closes until that customer 
enters service. ‘The gate then opens to admit all those who have 
accumulated outside the corral, and immediately closes again, re- 
maining closed until all those inside have acquired a server, and so on. 
Thus the customers are admitted in bunches to the corral, and once 
inside, they are served at random. It is assumed that the corral has 
an unlimited capacity. If there is no customer waiting when the gate 
opens, the gate remains open until there is again a customer who has 
to wait. 

A gated queuing system merits consideration, not only because of 
the anticipated improvement over random service, but also because 
it can overcome design problems in certain telephone equipment which 
might otherwise lead to poor service for some customers. For example, 
in some switching equipment, each server hunts in a fixed sequence 
over a group of customer-sources; when a customer is found, the 
server stops to provide service and then resumes hunting from that 
point. Such a hunting procedure may be desirable from a hardware 
viewpoint, but it can result in unequal service among customers. 
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Gates have been successfully used to improve service in the manner 
described above, that is, by temporarily blocking subsequent bids 
for service until all those customers present have been served. 

Some discussion of the gated queuing system is given in a 1953 
paper by Wilkinson.? Theoretical results were summarized in 1958 
by Beckwith,’ although he gives few details as to how the results are 
derived. The model we consider in this paper is more general than the 
one described above. We shall assume that as soon as there are k 
customers (rather than one) waiting, the gate closes until all k cus- 
tomers have entered service. Thus the parameter k can be thought of 
as a threshold at which the queue is judged to be too long to permit 
purely random service to continue; the system enters the “gating 
mode,’”’ and the gate admits customers in bunches as described pre- 
viously. If the gate opens and there is no customer waiting outside, 
then the system leaves the gating mode and the gate remains open 
until there is again a queue of & waiting customers. 

In Section II we obtain a recurrence relation for the probability- 
generating function of the nth bunch-size and, after that, the generat- 
ing function and moments of the equilibrium bunch-size distribution. 
Using these results we determine in Section III the Laplace-Stieltjes 
transform of the distribution of the equilibrium waiting time of an 
arbitrary customer. In Section IV we obtain the first two moments of 
the equilibrium waiting time, and we make some comparisons of the 
second moments of gated service for different values of the parameter 
k with those of random service and order-of-arrival service. 

In several places in the text, we specialize a general result to the 
case k = 1, since that is the simplest of our gated systems. This 
allows us to verify that our results then agree with Beckwith’s, and it 
often reduces a complicated expression to one that is more easily 
comprehended. 

We can assume, without loss of generality, that the system is 
initially empty. Since we assume that the traffic intensity, p = \/Nu, 
is less than unity, the number of customers will always return to zero 
in a finite length of time, regardless of how many customers may be 
present at the beginning. Consequently, our equilibrium results do 
not depend on the initial conditions. 

Several other authors have considered systems which operate in a 
manner similar to ours, but they all take k = 1 and assume that the 
underlying system is an M/G/1 queue. The “generations” in the 
M/G/1 queue, as described by Kendall‘ and later studied by Neuts,°* 
correspond to the bunches in our model. Nair and Neuts®:’ subse- 
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quently consider the waiting time distribution for the M/G/1 queue 
under the assumption that the queue discipline was either longest- 
processing-time-first or shortest-processing-time-first. 


II. DISTRIBUTION OF THE EQUILIBRIUM BUNCH-SIZE 


Let X, be the number of customers in the nth bunch; that is, there 
are N + X, customers in the system the instant after the gate closes 
for the nth time. Thus X,, is the number of customers waiting outside 
if that number is positive, and is k if there was no one waiting, when 
the gate was opened for the (n — 1)st time. Let 7, be the time it 
takes to serve X, customers. Then, denoting by p?(T, = t) the 
density of T,, at t, we can express the distribution of X, in terms of 
that of X,_1 by 


P(X, = j) = P(e =i) [P(X = j|Xea =i, Tra = 0) 
pU(T a1 = t|Xa-1 = t)dt. (1) 


But, given that X,»-1 = 7, T,-1 has a gamma distribution with param- 
eters 7 and Nu; and 


P(X, = j|Xa-1 = 14, Tri = 8) = P(X, = j|Traai = 6 


Oe 7 ee 
(iy ent + a, 4q i k. 


The extra term e—™ in this expression is the probability that no one is 
waiting when the gate opens; when this happens, the next bunch-size 
is necessarily k. Substituting these expressions into eq. (1) gives 


P(X.= i) = ¥ P&a=9(¢5) 


Cree) C75") +) 


where p = A/Nu and 6;, is the Kronecker delta. Now, letting 





P,(u) = v P(X. = j)ui 


be the probability-generating function of the distribution of X,, we 
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have 


PG) Pes (sa) h(i A Pe (<5): (2) 


Starting with Pi(u) = u*, we can determine successively the P,(u). 
When k = 1, eq. (2) agrees with Beckwith’s eq. (1). 

We wish to obtain the distribution of the equilibrium bunch-size 
X = lima. Xn. To see that an equilibrium distribution exists, we 
observe that the bunch-sizes X, form a Markov chain which is ir- 
reducible and aperiodic. Since we are assuming p < 1, the number of 
customers in the system cannot grow without bound, so the states are 
not all transient or null. Therefore all states are ergodic, and there is 
a unique equilibrium distribution.® The distribution of X = limn.. Xa 
has a probability-generating function P(u) = limi. P»(u), which 
satisfies 


This can be written in the form 
P[r(u)] — P(u) = Fu) 
by setting 
r(u) = 1/[1 + p(1 — u)] 
and 
F(u) = (1 — u*)PC1/(1 + 1) 1]. 


The solution to this functional equation is® 
P(u) = 9 — & Frau], 


where 1, is the nth iterate of r and n is a constant. To evaluate this 
expression, we first need to find r,(u). We have ro(u) = u and 


Sete 
1 + p — prn-1 


Tn (3) 
Letting y, be such that 
Th = Ynt — 1 + 1 
Yn p 
transforms eq. (3) into a linear homogeneous difference equation, 
which is easily solved. Thus we determine that 


1 — pu + (u — 1)p” 


tal) = Tou Fe — De 
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The solution to our functional equation is therefore 
so ates UR i eee =e |P 1 ) 
gC ae Et Feceacesraad ater 


Since P(1) = 1, we must have 7 = 1; since no bunch-size can be zero, 
P(0) = 0. This determines P[1/(1 + p) ], so that finally 


Ailterrto oll 





T— put (u— Ip 


P(u)=1-"°= 
Cc) 1 =: p” k 
n=0 {1 7 E i | 
1, hm 
h(0’ (4) 


where h(u) is the numerator in the second term of P(u) above. By 
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Fig. 1—Possible states of the system and the transition rates between states. 


0 IN SYSTEM 
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setting k = 1, we can reduce eq. (4) to Beckwith’s result, his eq. (2): 


Pu) = er | 





“—y 
p 
where 


p” 


1— a. 


gu) =u Y 


The moments of X are obtained by differentiating P(u). We find 
that 


si ace ae THO 
and 
wai ae a eh es |e 
Ve ©) = pgacylrestits maa | 


III. DISTRIBUTION OF THE EQUILIBRIUM WAITING TIME 


Let W be the waiting time to the point of entering service of an 
arbitrary customer when the system is in equilibrium. The distribution 
of W depends on which state the system is in when the customer 
arrives; i.e., it depends on the number of customers already in the 
system and on whether the gate is open or closed. These states, to- 
gether with the transition rates from one state to another, have been 
enumerated in Fig. 1. Let 


p; = P{j customers in the system ], J =20; 

pj = P[j customers in the system, gate closed _], 
7=N+1,N4+2,---,N+k —1, 

p; = P[j customers in the system, gate open |, 
g=N+1,N4+2,---,N+KA-1. 


Obviously, pj + p3 = p;. It is also clear that when j = N, the gate is 
open, and that when 7 2 N + k, the gate is closed. The values of the 
p; are just the equilibrium state probabilities of an M/M/N queue, 
which are 


A/p)? N) , 
= Ny = on. Po, fas Oy A Sot gv 


; .(pN)*% 
PNn+i = PPNn+j-1 = p’pn = p’ ew) Po; j>0 


_ Pf Xe! (pN)s (oN) 7 
pom | jt Wid | 
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To find pj, we equate the rates at which the system leaves and enters 
the state {7 customers in the system, gate open}. From Fig. 1, we see 
that 


Chee P)PN-+43 = ppy+j-1 ae PN-+i+19 j=H1,2,°°+,k— 1, 


with pyiz = 0 and pi = py, where py is known from the above. 
The solution to this equation is 


J — pk 


Dhets = Tae Pm j3=0,1,°--,k-1. 
It then follows that 


Pu+i = Pn+i — Duss 


1 — p’ : 
= Tae tPw, gul,2,-+-,k-1. 
To find the distribution of the waiting time W, we shall consider 
what happens when an arriving customer encounters one of the 
following conditions: 


Hi: < N customers in the system, 
H: the gate is closed, 
H3.;: N + 7 customers in the system, and the gate is open, 





7=0,1,---,k-1. 
From the above computations, 
eo ; 1 
P(H1) =1— 2% peppy = 1 — > byw, (5) 
j=0 —p 
k-1 eo kp* 
P(H2) = DO prag+ Do pws = T— 9 P™ (6) 
j=l jak p 
and 
pi — p* ; 
P(H3.;) = PN, 7 = 9, 1, 98 glee (7) 


1 — p* 


An arriving customer who encounters H; immediately gains access 
to a server, So 
1, t=0 


0, %$t>0 (8) 


P(W = t|H) -{ 
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— ——-G,+- ————— 
Pe eee Gene eet 
LAST CUSTOMER GATE CUSTOMER 
GATE-CLOSING ARRIVES OPENS ENTERS SERVICE 


to 


Fig. 2—Relations among the variables used. 


An arriving customer who encounters H», will have to spend some 
time, Z, outside the gate waiting for the gate to open (see Fig. 2); 
once he gets inside, he will then have to spend some additional time, 
Y, waiting to be chosen for service from the bunch he entered with. 
The total amount of time the customer spends waiting to enter service 
is therefore W = Z + Y, where Z can be regarded as the residual 
lifetime of an interval G» during which the gate is closed. Thus we 
can write 


p(W = t|H») is ” p(Z+Y =t, Go = y, Z = a)dyde 
az=0 /y=0 


[Of wer =6- 2lG= 9,2 =2) 
=0 y=0 
-pi(Z = x|Go = y)p*(Go = y)dydx. (9) 


We first need the distribution of Go. Disregarding our arbitrary cus- 
tomer temporarily, let G be the equilibrium length of the time in- 
terval during which the gate is closed. If we are given that 7 customers 
entered the waiting room when the gate last opened, then the gate 
will remain closed for j service times. Since the equilibrium probability 
that 7 customers entered when the gate last opened is P(X = 9), 
we have 


pi(G = y) = 2 P(X = jb*(y), 
= 
where b(y) = Nu exp (—Nuy) and the asterisk denotes convolution. 
The mean of this distribution is easily found to be 


E(X) k 


BG) Nae = DAO) 


Let to be the epoch at which the arbitrary customer arrives, and let 
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Go be the length of the G-interval containing t). Then (see Ref. 10) 
p4(Go = y) = aa PMG = 9) 
E(G) 


1 e oe 
= 7 Nay(1 — p)h(0) u PX = NoMy). (10) 
= 
Next, the distribution of Z, given Gp = y, is uniform on (0, y), so 


7 S44 
pi(Z = x|G = y) = 49 (11) 
0, 


u> Y. 
We also have 


p(Y =t—2|Go=y, Z = x) = p'(Y =t — 2|G@o = y) 


8 


= > p4(Y = t — x|n arrivals in (0, y), Go = y) 
1 


n 


-p?(n — 1 other arrivals in (0, y)|Go = y) 


ar (= beg — »)) . oS ene, (12) 


where the first factor in each summand reflects the fact that the ar- 
bitrarily chosen customer has probability 1/n of waiting one service 
time, 1/n of waiting two service times, --- , and 1/n of waiting n 
service times. Substituting eqs. (10), (11), and (12) into eq. (9), we 
obtain 


p*(W = t| He) 
= Nal —o)hO) EE PK =| fl we - a) 
jg=in=1 l=1 J/2=0 
fo pti(y)e re OO dyda. (13) 


Notice that the range of integration of the inner integral was reduced 
from (0, ©) to (z, «) because of eq. (11), and that of the outer in- 
tegral was reduced from (0, ~) to (0, £) because b*'(é — x) = 0 for 
x >t. By setting k = 1, Nu = 1, and d = p in eq. (13), we can ob- 
tain Beckwith’s expression for the corresponding density in his model, 
the density of W given that the arbitrary customer finds more than 
N customers in the system. 

We can calculate the Laplace-Stieltjes transform, ¥(s), of the dis- 
tribution (13), but the algebra is long and tedious. The results are 
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given below, one in terms of P(X = 7) and the other explicitly: 
v (Nus) 
_ (1 = p)h(0) P(x = 1) log Lt 9 + 8 + 80) 





pks? (l+s+p)(1+s)— »p 


2 | = 1 \4 7 its \7 
+ BaP = al1- (7p) Gait.) 


as (a aa eae s) — 3) || (14) 
ree ke — | (Si Ne 


(1 +s — p — sp™)(1 +8 — p — sp™) 
-log ( (1 —p)[A+s)? —~p—s(l+s+ =) 


ECCLES) (nite) 
ae 


(L— p)E(+s)?+ sp] \" 
a eercriceesrc El 




















= 1 Aah (1 + s — p — sp*)* 

1 — ptt p (Il +s —p— sphyei 
_ il _ ee Seep )e 

1l+s (+s — p — sp**)F1 

Le Eas! Sp Ss ek pel |: (15) 
Lees: LO 6)t = pel scr pp ir 
When & is set equal to 1 in eq. (15), the sum from 2 to k is trivially 
zero, and the expression in the last pair of braces collapses to zero; 
the transform then reduces to 
Chip) p” 

ps* nso (1 — p®*)? 
X lo (leap Sep) Te pee) 

PL @— eG +s? — 0 — sd +s + aor] |’ 
(k = 1). (16) 








+ 








¥(Nus) = 
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The final portion of the waiting time distribution needed is 
f(t) = p?(W = t|H3.;), 7 = 0,1, ---, k — 1. Let W; be the waiting 
time of a customer who arrives to find the gate open and N + j in the 
system. Then W; has the density f;(¢), and, letting H,(t) denote the 
exponential density \e~™, 


f, = X N+ Nu Eyanu* fi-1 + —_ x Mi 


| Bae + aren sage 
= 0, ,k—-2 


1 1 1 
fra = j exw + 7 Ene + nas + 7 EM 


The reasoning behind these equations is that, if an arriving customer 
finds the gate open and 7(< k — 1) customers waiting, then either 
the next event is an arrival, in which case he waits until that event 
plus an additional time distributed as W;,1, or else the next event is a 
departure, in which case he waits until that event, after which either 
he is served immediately [with probability 1/(j + 1)] or someone 
else is chosen and our customer must wait an additional time dis- 
tributed as W;_1 [with probability 7/(j7 + 1)]. If, on the other hand, 
the arriving customer finds the gate open and k — 1 customers waiting, 
then the gate shuts behind him and he waits either 1, 2, --- , or k 
service times, each having probability 1/k. Taking Laplace-Stieltjes 
transforms of the equations, and denoting the transform of f,(¢) by 
¢;(s), we obtain 


oj¢;(Nus) = (1 + p + 8)76;-1(Nus) — (7 — 1)¢;-2(Nus) — 1, 


p= 1,-+-,k—-1, 
oxa(Nus) = ae , 4 (a7) 


For any particular k, this set of equations can be solved explicitly by 
successive substitution. Finally, using eqs. (5), (6), (7), and (8), we 
can represent the Laplace-Stieltjes transform of the distribution of 
the waiting time W as 





6(Nus) = 1 — p> — pw + pw PE os(N us) 


a8 pw¥(Nus), (18) 





kp* 
1 — p* 
where ¢;(Nus) is given by eq. (17) and ¥(Nus) is given by eq. (14) 
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or eq. (15). When k = 1, eq. (18), with the help of eq. (16), can be 
written explicitly as 


1 (1 — p)? = p” 
Las punt 2 PN o> Gd — p=) shy 


(1 + s — p — sp"*)(1 + 8 — p — sp)**? = 
oe Goto |) =D 








1 
$(Nus) = 1 — reset Jaa 





IV. MOMENTS OF THE EQUILIBRIUM WAITING TIME 

Since the mean waiting time does not depend on the queue dis- 
cipline (see Ref. 11), the mean is the same as for a simple queue with 
service in order of arrival, i.e., 


_< os eae ee Py 
ene DyPre’ Nu Nu(l — p/? 





The second-moment computations are fairly lengthy, finally yielding 
2pNn 
E(W?) = ——>——~ [M'(p) + M(p) — p*3M'(1) — pF OM (1 
(W?) Madey (p) (o) — p (1) — p (1) J 


kpnp*(1 — p) | (k — 1)(k — 2) |? + 2p + | 
(Nu)?(1 — p*) 6 t= p* 


ae eS 19) 
(I. = p?)(L —"p*) (i= p) Gp). = p2). | 


where M is a function defined by 


+ 


+@-D| 


MO = 6; (0) = Foe W | Hs.) 


The variance of W is obtained by subtracting the square of the mean 
of W: 


Var (W) = E(W?) — Re: = 


For comparison purposes, we also need the second moments for 


order-of-arrival service and for random service. When service is in 
the order of arrival, we obtain" 


E(W?) = E(W2|W > 0)P(W > 0) = Wore ys (20) 


When service is at random, the second moment can be written as" 


ee ee 
(Nu) (l = p)* 2° p 





E(W?) = (21) 
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Observe that the second moments depend on N only through the 
factor py/(Nu)*, assuming the value of p is fixed. This is true also 
for the second moment in eq. (19), since each of the M-terms con- 
tains a factor (Nu)—!. It is therefore convenient to consider the ratio 
of (W?) for the gated system [eq. (19) ] to the second moment for 
the order-of-arrival system [eq. (20) ], i.e., 


R(ky p) = Net ee [pM"(p) + M(p) — pkM"(1) — p*M(1)] 





kor (=p) | (k= 1) = 2) k + 2p + 2] 
2(1 — p*) 6 =p? 


“a ae | 2+ 4p + 5p? + 2p? + pt 








(1 =p?) (=p?) (l= pil 97)(1 —‘p*) 


Because this ratio is independent of N, it provides a useful tool for 
examining the effect of the value of k on the second moment of the 
waiting time W. Thus we shall be interested in determining its 
properties. 

The function R(k, p) has been plotted in Fig. 3 on the interval 
0 <p <1 for a variety of values of k. We observe that R(k, p) is 
bounded from below by unity, increases as k increases, and is bounded 
from above by 2/(2 — p), the ratio of eq. (21) and eq. (20), which 
corresponds to k = . This general behavior is, of course, just what 
we expected a prior7. In order to demonstrate analytically that 


+@-1| ri (22) 





1=R(k, ep) = (23) 


2—p’ 
we first introduce the inequality 

Af S 1j+2 2 

Nu 2 Nu 2 2-p’ 
j3=0,1,---,k-—1. (24) 





= E(W|H3.;) S 


The left half of this inequality is demonstrated by considering what 
happens when an arriving customer finds 7 others waiting, but no 
more arrivals are permitted to enter the system: the customer’s 
expected waiting time will decrease to (7 + 2)/2Nu, since the customer 
will have to wait either 1, 2, --- , or 7 + 1 service times, each with 
probability (7 + 1)—!. The right half of the inequality is demonstrated 
by considering what happens when there is no gate at all to block 
future arrivals when a threshold k is reached: the mean waiting time 
Ei (W|#3.;) would then increase to the corresponding mean in a system 
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R (k, p) 





Fig. 3—Ratio of second moments. 


employing purely random service. In a random service system, a 
customer who arrives to find j others waiting has an expected delay of 
ee 


Nu 2 2-—p’ 





PSO Ae he A 


a fact which is derived in the appendix. 
Using the left half of eq. (24), we can substitute (7 + 2)/2Nu for 
E(W |#H:3.;) in R(k, p) and combine terms, obtaining 
R(k, p) 21 
kp**(1 — p) | 3k(1 — p)(1 — p?*) k?(1 — p)? 
12(1 — p*) (E+ eee): Trae 
from which it is obvious that R(k, p) 2 1. Similarly, using the right 


half of eq. (24), we can substitute (7 + 2)/Nu(2 — p) in Rk, p). 
Combining terms, we obtain 


R(k pe gs ee 
rh =2—p @-e)G—r) | 120 +0 + 0%) 


k(1 — p)(2 + 2p + 5p? + pt) , 2+ 7p + 10p? + 8p* + 3p! 


BRE py (hp ee) a eae V6 eT 
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from which it is clear that R(k, p) = 2/(2 — p). Thus, eq. (23) is 
established. 

Perhaps the most striking feature of Fig. 3 is that all the curves 
approach the same value, 7/6, as p— 1. It is easy to demonstrate, 
by using eq. (22), that this must occur. Since the conditional means 
E(W|4H3.;) remain bounded as p—1, M(p) — M(1) and both are 
finite. Thus the first term of R(k, p) approaches zero as p — 1. In the 
second term, the factor (1 — p)4 in front causes all but the last term 
in braces to approach zero. Thus, 
kp} 2+ 4p + 5p? + 2p? + p* 7 





as repel tpt» +p) “1+ p)iteote) 6 

In order to gain some insight as to why the curves meet at a common 
point at p = 1, we will find it helpful to consider a supplementary 
variable, the fraction of time the system spends in the gating mode. 
When the system is in equilibrium, this is simply the probability that 
an arriving customer finds the gate closed, and is given by eq. (6): 


kp* 
Lp 





P (system is in gating mode) = 7 PN: 

This quantity has been plotted in Fig. 4 as a function of p, for the 
arbitrarily chosen value N = 7. It can easily be seen (and is intuitively 
obvious) that when p is very close to 1, the system spends almost all 
its time in the gating mode. But when the system is in the gating mode, 
the system’s operation is independent of the value of k; it is only 
when the gate is open that the threshold value k can have any effect. 
Thus, as p approaches 1, the system becomes independent of k; so we 
can expect the curves in Fig. 3 to be independent of k at p = 1. 

Another feature of the curves in Fig. 3 is that the slope at zero for 
k = 2 is the same as the slope of 2/(2 — p); but the slope for k = 1 
is zero. The reason for this becomes clear when we realize that when 
k = 1, order-of-arrival service is guaranteed until there are N + 3 
customers in the system, while k 2 2 makes it possible for passing to 
occur as soon as there are N + 2 customers in the system. For small 
p, N + 3 in the system is much less likely than N + 2. 

The ratio in eq. (22) is convenient because it is independent of N. 
The variance of a distribution, however, is also frequently of interest. 
It is clear from Fig. 3 that the second moments, and therefore the 
variances, increase with increasing k. There is, therefore, a similar 
family of curves, one family for each value of N, obtained by computing 
the ratio of the variance with threshold k to the variance with order-of- 
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Fig. 4—Equilibrium probability that gate is closed (N = 7). 


arrival service. In Fig. 5 we have plotted some of these curves, again 
for N = 7. It is obvious that these curves have the same general 
shape as those in Fig. 3. The main differences are (2) that the ratio of 
the variances when k = © (random service) goes to 3 as p — 1, while 
the ratio of the second moments when k = © goes to 2, and (72) that 
the ratios of the variances for k < © go to 8/6 as p— 1, while the 
ratios of the second moments go to 7/6. These facts are easily verified 
analytically by letting p — 1 in the actual expressions for the ratios. 

Since the second moments (and variances) increase as a function of 
k, there arises the question of why one might consider a threshold 
value k greater than 1. Clearly, if the queuing systems were otherwise 
equivalent, one would prefer k = 1 to any larger value of k. But it is 
equally possible that a queuing system would be more costly to operate 
when it is in the gating mode, since more bookkeeping is necessary : 
each waiting customer must be classificd as “inside” or ‘‘outside”’ 
the waiting room, and these labels must be changed when the gate 
operates. One would then prefer a system in which the gate were used 
as little as possible (large k), consistent with an acceptable quality 
of service. The resulting tradeoff between cost and quality of service 
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2.2 


2.0 





Fig. 5—Ratio of variances (N = 7). 


can be resolved only by examination of the particular application at 
hand. 


V. CONCLUDING REMARKS 


In summary, our main result is the specification of the distribution 
of the equilibrium waiting time of an arbitrary customer in a queuing 
system whose service discipline is a compromise between service in 
order of arrival and random service. Our model contains a parameter, 
k, which determines how ‘‘close’”’ the discipline is to order-of-arrival 
service or to random service. We have seen that the variance of the 
waiting time is bounded from below by the variance for order-of- 
arrival service, and that as k increases, the variance increases, ap- 
proaching the variance of the waiting time when random service is 
employed. We also found a convenient quantity, R(k, p), which is 
independent of the number of servers, and which, together with Fig. 3, 
allowed us to examine the effect of the threshold & on the waiting time. 
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Figure 3 shows how, for fixed k > 1, the service changes from ‘“‘nearly 
random” to ‘‘not quite order-of-arrival”’ with increasing load, and how 
this transition occurs.at higher loads as k increases. 

There are, of course, other variables which can be derived from the 
system we have described; for example, in studies of equipment life, 
it might be useful to know the distribution of the number of gate- 
closings that occur in an interval of length ¢. It did not seem worth- 
while to pursue such questions in the present study, which deals with 
gating from the viewpoint of traffic performance. 


APPENDIX 


Suppose we have an M/M/N queuing system employing random 
service, with arrival rate \ and mean holding time y~!. The mean wait- 
ing time to the point of entering service of an arbitrary customer in 
such a system is 
___ Py 

Nu(l — p)?’ 
where p = A/Nu <1 and py is the known equilibrium probability 
that there are exactly N customers in the system. We wish to deter- 
mine the mean waiting time, m;, of a customer who arrives to find 
N + j other customers in the system. The m,’s satisfy 


= 1 ON Nu i] 
M+ Nu \+Nu SN Gg sed 


The rationale for this equation is that a customer must wait at least 
until the next change of state; the mean of this initial delay is 
(\ + Nu). If the change of state is caused by an arrival [which 
occurs with probability \/(\ + Ny) ], then the customer will have to 
wait an additional period of time whose mean is m,;+1. If, on the other 
hand, the change of state is caused by a departure [ which occurs with 
probability Nu/(\ + Nu) ], then with probability j/(j +1) our 
customer will not be chosen from the group of 7 + 1 customers, and 
he will have to wait an additional period of time whose mean is mj-1. 
We now introduce the function 


m 


~ 


mM; ~+ Miri + M1, j 20. (25) 


H(z) = myer. 
= 


This series converges for x = p, since the mean waiting time of an 
arbitrary customer is 


m= pnH (p) = Nall Sp? (26) 
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Multiplying eq. (25) by (1 + p)(j + 1)a’ and summing on Jj, we can 
obtain a first-order differential equation: 


L+p—2 all 1 . 
t=9G—) °° N= eG =e) 


The solution to this equation is 


1 = —=\" 2 —7 
p Nu(1 — x)?(2 — p) 


H’(a) + 





H(x) = C(1 — x) € 
Using eq. (26) as the boundary condition, we see that C' must be zero 
in order that H(x) remain finite as « — p. Thus 


Soke 8 as 
Null =—aZ)2 =p) 


The power series expansion of this function is found to be 


H(x) = 


-) + 2 
Hie) = i] 
(*) » Nu(2 — p)” 
therefore, the means we desire are given by 
pi TES 
™ ~ Nu(2 — p) 
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Losses and Impulse Response of a 
Parabolic Index Fiber With 
Random Bends 


By D. MARCUSE 
(Manuscript received March 29, 1973) 


The coupling coefficients of the modes of a parabolic index fiber with 
randomly curved axis are derived and are used to compute its excess losses 
and impulse response. It is found that bends with a period comparable 
to the natural ray oscillation period in the parabolic index medium are 
catastrophic. The average radius of curvature R, of a guide composed of 
circular sections with an average length of 1 cm must not decrease below 
approximately R, = 1 m. Mode coupling by random bends has the tend- 
ency to reduce the width of the impulse response function. However, this 
improvement ts accompanied by losses. Reducing the width of the ampulse 
response for coupled mode operation to half its uncoupled width causes 
0.7 dB additional loss, a ten-fold reduction of the pulse width costs 18 dB. 


I. INTRODUCTION 


Optical fibers with parabolic refractive index profile’? have less 
pulse delay distortion than conventional fibers with piecewise con- 
stant, discontinuous index distribution.? The width of the impulse 
response increases in direct proportion to the length of the fiber. It 
is well known that an improvement of the impulse response results if 
the modes are coupled among each other.*:> In the presence of mode 
coupling the width of the impulse response increases only propor- 
tionally to the square root of the length of the fiber. 

Pulse propagation in multimode parabolic index fibers is studied in 
this paper by means of converting the coupled power equations to a 
partial differential equation.*-’ Random changes of the direction of 
the waveguide axis are considered as the coupling mechanism. 

The problem is simplified by assuming that the modes of the para- 
bolic index fiber are essentially the same as the modes of an infinitely 
extended square-law medium. The fiber boundary is included in the 
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description by requiring that modes interacting with it suffer high 
losses so that an effective cutoff exists. Modes below the cutoff value 
propagate as if they were in an infinitely extended medium. At cutoff we 
demand that the modes do not carry power. The effect of the wave- 
guide wall is thus taken into account as a boundary condition that 
has to be satisfied by the solutions of the partial differential equation. 

This formalism provides information about the width of the impulse 
response function and the losses associated with the coupling mecha- 
nism. The achievable improvement in pulse width due to mode coupling 
can thus be expressed in terms of the associated loss penalty. We 
conclude that mode coupling is capable of improving the already 
favorable impulse response of the parabolic index fiber. However, 
this improvement of the width of the impulse response is accompanied 
by excess losses. The product of the square of the pulse width ratio 
(width of the impulse response of coupled modes to the uncoupled 
pulse width) times the loss penalty is independent of the waveguide 
parameters and the statistics of the axis deformation. There is thus 
no hope of reducing the loss penalty of delay distortion improvement 
by optimizing the waveguide parameters. 


II. MODES AND COUPLING COEFFICIENTS 


We use the modes of the infinite square-law medium. The refractive 
index distribution is assumed to be of the form 


nt = ng (1-285) (1) 


The fiber radius is at r = a. However, the modes are assumed to be 
unaffected by the fiber boundary if their mode number remains below 
a cutoff value. We use linearly polarized modes and obtain for the y 
component of the electric field® 


(fr) m2) a) om 
€0 WwW Ww = 
Bg = tBz, 


(nor2?t 2p! g!)w 


(2) 


There is also an electric field component in axial direction. But, for 
small refractive index changes, it is negligible. The functions H, and 
H, are Hermite polynomials of order p and q, the radius r is defined by 





r= e+ y, (3) 
and the mode radius w is 
v2a \3 
2 (sae ) 4) 


MULTIMODE PARABOLIC INDEX FIBERS 1425 


with the free-space propagation constant 


k = wVeouo. (5) 
The propagation constant of the mode is given by the expression® 
2vV2A i 
g=nk]1—-otaty]: (6) 


The orthogonality of the modes and their normalization follows from 
the following equation: 


ae [ i Ep gH yydxdy = PSpp'Sqq'- (7) 


The asterisk indicates complex conjugation. 
The cutoff condition for the guided mode has been derived in Ref. 9. 
The permitted maximum values of p and q are defined by the relation 


A 2 
(p+qe= fh naka = — (8) 


Any deviation of the parabolic index fiber from its perfect geometry 
can be expressed by a change of its refractive index distribution. 

Changes in the direction of the waveguide axis can be expressed by 
the following index distribution: 


nt =n {1 ~ 23 LGe — faye + vl} () 


We consider waveguides bends in only one plane for simplicity. The 
results thus obtained can easily be extended to the general case. The 
appropriate coupling coefficient for this type of index change is 

K pq,p'¢! =, AP (Bag = Bora) oe _ dz E508’ a dxdy. (10) 
Bends of the waveguide axis couple even modes to their immediate 
odd neighbors and odd modes to their immediate even neighbors. 
Only the following coupling coefficients are different from zero: 








df 
kwA dz .nokwd 
Koq.pttia = — ND pe gia al — Vp f(2). (11) 


Because we restricted the waveguide curvature to the x — z plane 
there is no coupling between the modes with different values of g. 
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All coupling coefficients with different gq values vanish. The derivative 

of f(z) was replaced by the function itself with the help of the relation 
‘a: 5 ee 

This replacement is permissible since it is the spatial frequency 6 — 8’ 

of f(z) that is responsible for the coupling process. 


III. COUPLED POWER EQUATIONS 


Pulse propagation in multimode optical fibers can be described by 
the following set of coupled equations for the average power P, 
carried by the modes :® 


oP, , 1 OP, 


0z Vy, Ot 








N 
=~ a4Py + Y hw(P, — Py). (13) 


The single label v indicates both p and qg. The power coupling coeffici- 
ents h,, are defined by * 


hus = Roy = |Kuo|? F(6, — B,). (14) 


The coupling coefficient (11) enters the power coupling coefficients 
via the definition 
Ky = Ky f(2). (15) 


The power spectrum F'(6, — 6,) is defined by the equation 


F(6) = ( = i ” t(e)e-ttede ». (16) 


The symbol ( ) indicates an ensemble average. 

Coupling between the modes of the parabolic index fiber is an ideal 
application for a diffusion theory of the power coupling process.’ 
Since only nearest neighbors couple directly to each other, power 
redistributes itself by jumping from mode to mode in the same way 
as particles diffuse through real space. If the mode number is very 
large we can consider the set of discrete modes as a quasi-continuum 
and change the equation system (13) into a partial differential equation. 
To accomplish this transformation we consider the following expres- 
sion using hy, = Avy: 





N 
2X Ry(P, _ P,) = Rust (P pst = P,) = hurl Pu a P24). (17) 
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Considering yu as a continuous variable » = 6 we use the approximation 


hutisu (Pest ~ Pi = hei a =: Py-1) 


= As {H(o +00) (Se) - ho (3), 


= (ao 2 (1 =): (18) 


With Aé = 1 we can write (13) as a partial differential equation 


ap , 1 aP 
Vv 


= a /, oP 
= at =-aP+5(h }. (19) 


00 
The propagation constant (6) can be approximated as follows: 


p= nk - Wt g+y (20) 


with p = 6. The difference of the propagation constants of adjacent 
modes, 


AB = B(0 + A) — B() = — “““ = -Q, (21) 


is independent of 6. The power spectrum entering (14) contributes to 
the coupling process only at one fixed spatial frequency and is in- 
dependent of the variable 6. The power coupling coefficient (14) can 
be expressed with the help of (4), (11), and (21) as follows (p = 86): 


v2n ok A? 
ae 


h(0) = F(Q)6. (22) 





With (22) the partial differential equation (19) assumes the form 


aP | 16P _ 


+57 TP +K| (23) 


oz | ov at a6 86 


2 
9 eP . oP | 
with 
_ V2nokA} F(Q). 


K a 


(24) 


We assume that the attenuation coefficient a in (19) is constant and 
describe the high loss, that we must attribute to modes interacting 
with the waveguide boundary, by means of the boundary condition 


P(z, t, €) =0 for 6 = 6c. (25) 
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The cutoff value 6 = 6, follows from (8): 


9 _@e? _ nokavA _ 





The slope 0P/06 determines the rate of power diffusion. Since no 
power can be lost at 6 = 0 we must also require 


oP 


IV. STEADY-STATE POWER LOSS 


We begin the discussion of the solutions of (23) by neglecting the 
fact that each guided mode has a slightly different group velocity and 
consider v(@) = const. We construct a solution of (23) by introducing 
the trial solution 


P(z, t, 0) = e~+%G (6). (28) 

Substitution of (28) into (23) yields the ordinary differential equation 
OG dE 6 i pe: 

67@ tag t Re = 2. (29) 


The normalized solutions of this equation that satisfy the boundary 


conditions (25) and (27) are 
J , ac ) 
1 ° if 6. (380) 





G, 6 = 
re ahah) 
with 
u 7 Oe. (31) 


The parameters u, are determined as the roots of the equation 
Jo(uy) = 0. (32) 
The functions G,(@) are mutually orthogonal. 


i * ,()G,(6)d0 = Sy. (33) 


The general solution of the power equation (23) is obtained as the 
superposition of the trial solutions 


P(a, t, 0) = e- > ¢,G,(0)e7™™, (34) 
v=1 


The expansion coefficient c, can be determined from the power dis- 


MULTIMODE PARABOLIC INDEX FIBERS 1429 


tribution at z = 0 with the help of the orthogonality condition (383), 
Bc 
ei = 7 G,(6)P (0, t, 6)d0. (35) 
0 


The eigenvalues oc, are obtained from (81) and (82), 


2 
_ Kuve 


an ae. (36) 





The eigenvalues increase with the increasing values of the roots uy. 
It is thus apparent that only the first term in the series (34) needs to 
be considered for large values of z. The steady-state power distribu- 
tion is thus described by the equation 


Pz, t, 0) = cre t9*G, (0) for 27, (37) 


After an initial transient has decayed, the power distribution (versus 
mode number 6) assumes the steady state (37). The power loss in the 
steady state is the sum of the constant loss a, that was assumed to be 
the same for every mode, plus the loss value o; that stems from mode 
coupling due to waveguide curvature. With 


uy = 2.405 (88) 
we obtain the steady-state curvature loss from (24), (26), and (86), 


i 2.045 + rok A? (9) 
. nokavaA _ ) P (39) 
ia) A 





Because of our assumption that the waveguide is curved only in one 
plane the steady state losses depend on the mode number gq. For small 
values of gq we find low curvature loss 

_ 2.894 


oO, = 
a} 





F(Q) for q=0. (40) 


With increasing values of g the losses increase until they reach in- 
finitely high values. 

In any actual cases it is unrealistic to assume that the waveguide 
would be bent in only one plane. Bends in the perpendicular plane 
couple modes with different gq values. The total steady-state loss is 
thus a weighted average of the losses (39). The weight factor is the 
number of modes for each value of g. According to (8) we have 


p= N(@Q) = 2 naka —q (41) 
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different modes for each value of g. The average loss that results from 

coupling all the modes by random bends in both planes is thus (see 

appendix) 

5.84 
a‘ 





4 nokaVA/2 
A-cion f  wN@a= FO). — (42) 

Comparison with the loss coefficient (40) for the mode group with 
q = 0 shows that the total fiber loss (42) is just twice as large. The 
loss coefficient (40) is representative of a slab waveguide. The actual 
loss of the round fiber can thus be deduced from the slab waveguide 
model. That the loss coefficient of the fiber is twice as large as the slab 
waveguide loss might be expected, since F'(Q2) stands for the ampli- 
tude of the power spectrum for bends in only one plane. The fiber is 
assumed to be bent in both planes with equal power spectra. The 
effect of both bends add, doubling the loss coefficient. 

In a fiber cable the fibers may (or may not) be twisted around each 
other. Such twists could introduce an almost sinusoidal deformation of 
the fiber axis. The length A of the period of sinusoidal deformations 
that should be avoided follows from (21), A = 27/Q. For numerical 
estimates we are using a fiber with radius a = 4.85 X 10-3 cm and 
A = 0.014. The critical period for this fiber is A = 0.18 em. If such a 
period should be built into the fiber by the method of cable construc- 
tion we can estimate the losses that would be caused by a given ampli- 
tude. F(Q) has the dimension of cm’. It can be interpreted as the ratio 
of the square of the amplitude of the sinusoidal deformation and the 
spatial bandwidth that may be caused by random phase changes. 
With Q = 34.5 cm! let us assume a spatial bandwidth of AQ = 3 
cm, An excess loss of é; = 10 dB/km = 2.3 X 107-5 cm™ would re- 
quire F(Q) = 1.57 X 10-® cm’. The square of the amplitude A is 
given by the product of F(Q@) with the spatial bandwidth. We thus 
obtain the amplitude of the sinusoidal deformation of the fiber axis 
that causes an excess loss of 10 dB/km: A = [F(Q)AQ]}} = 6.9 
X 10-7 cm = 69 A. Sinusoidal axis deformations at the critical wave- 
length are thus seen to be extremely dangerous. 


V. LOSSES FOR A STATISTICAL MODEL 


Even though it is only one spatial frequency of the power spectrum 
F(6, — 8.) that determines the steady state loss (42), it is hard to 
guess the amplitude F(Q) that might be expected at this spatial fre- 
quency Q. In order to gain insight into the expected steady-state 
curvature losses it is necessary to consider statistical models. We 
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consider a model consisting of waveguide sections with constant 
curvature whose magnitude and sign varies randomly. 
The power spectrum can be written as follows: 


F(Q) = 7 «| if ” F(z)e-iMede » 
= ar ( i PY ins ). (43) 


dz? 

The step from the function f(z) to its second derivative involved two 
partial integrations. The end points of the integration range do not 
contribute if we assume that the randomly disturbed guide is con- 
nected to two perfectly straight waveguide sections so that f(z) and its 
first two derivatives vanish at z = 0 and at z = L. 

For waveguides that are only slightly bent we can consider the second 
derivative of f(z) as the curvature function 1/R,(z). We denote with 
C(u) the autocorrelation function of the curvature function, 











C(u) = (44) 


(ROR ETD) 


It is well known that the power spectrum of a function is equal to the 
Fourier transform of its autocorrelation function.!2 We may thus write 


1 ; 
F(Q) = 5 ie C(uje*@4du, (45) 
The autocorrelation function of waveguide sections with piecewise 
constant curvature and fixed length D is 
D = | | 2 


GG) = D juj sD 
0 |u| > D. 


(46) 


The parameter x? is the variance (square of the rms value) of the curva- 
ture 1/R.. Substitution of (46) into (45) results in 


F= 





2K? 

Doe (1 — cos QD). (47) 
An average over this expression, that allows us to consider D as an 
averaged quantity, leaves us with 

2k? 


F(Q) = pes 


(48) 
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From (21), (42), and (48) we obtain the following expression for the 
steady-state loss of our statistical waveguide model: 


272 
sare 1.474; (49) 


D2? 





As a numerical example, we consider a parabolic index waveguide 
with the following parameters: 


\ = 1 um, free-space wavelength (k = 6.28 X 104 em7?) 

a = 4.85 X 10-3 cm, waveguide radius 

No = 1.56 

A = 0.014 . (50) 
w = 7.69 X 10-4 cm, mode radius 

ee 

w 


With these data we have (x in cm~, D in cm) 
K2 
oY) = 0.17 D em}. (51) 


We may now ask for the rms value of the curvature that is required 
to cause a steady-state loss of 2.3 X 10-> cm = 10 dB/km with an 
average length of the waveguide sections of D = 1 cm. We find from 
(51) « = 0.0116 cm~ or 1/x = 86 cm. Our result tells us that a wave- 
guide composed of individual sections of constant curvature of average 
length D = 1 em with R,= 1m radius of curvature has 10 dB/km 
additional loss. For our derivation we assumed that the high-order 
modes suffer very high losses since their fields reach into the vicinity 
of the waveguide wall. Whether the interaction with the outer wave- 
guide boundary causes high losses depends on the construction of the 
waveguide. If the outer surface is rough or coated with an absorbing 
material to reduce crosstalk, the losses are high and our estimate 
applies. 


VI. PULSE DELAY DISTORTION 


It was shown in Ref. 5 that the width of the impulse response of a 
multimode fiber with coupled modes is given by the equation 


At = 4VpL. (52) 


Lis the length of the waveguide and p is the second-order perturbation 
of the eigenvalue o; defined by (36). For the discrete case we write 
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G,(6) = G,™ and obtain p in the form 


N 1 1 2 
| oe (; = ~) CG, | 
p=1 \Up Va . (53) 


2 Oy — O1 


Mz 


p= 


vy 


The average group velocity v, actually does not contribute to (53) 
on account of the orthogonality of the vectors G™. With the assump- 
tion of a continuum of modes we obtain instead of (53) 


oe = | ln ge (ae : ~) G(0)G,(0)a8 | ; (54) 


The inverse group velocity is obtained by approximating (6), 


25 i 
BR nk — —— (p + g) — TE (p+ 9, (55) 
and taking the derivative. With v, = c/n, we obtain (with p = 6 and 
c = light velocity in vacuum) 

1 1 1dpn A 


v(6) Vo c dk Cc cn,k?a? (6 + 26g ote q’). (56) 


The term with q? does not contribute to the following integral because 
of the orthogonality relation (33): 


[, (0 = ~) G1 (6)G, (0) d8 
~ ant Ea - Td | +a}: (57) 


enok?a? (u;, — ui)? (u; — ui)? 


Each mode group with a given value of q has a different spread of 
the group velocities of its uncoupled modes. However, we have seen 
that the waveguide losses could be obtained from the simpler slab 
waveguide model. This simplification is expected to apply also to the 
pulse distortion problem. The slab model is obtained by setting g = 0. 
The spread of inverse group velocities is largest for g = 0 since the 
allowed 6 range is largest in this case. However, even though the spread 
of the group velocities is reduced for increasing values of qg, the mode 
groups with different q values arrive at different times. This delay 
distortion is reduced by coupling of the different mode groups by 
means of waveguide bends in the perpendicular plane (perpendicular 
to the plane coupling the modes with different p values). The mode 
group with g = 0 can be excited by shining light into the fiber that 
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is collimated in one plane but spreads in the plane of the bends in 
such a way that all modes with g = 0 and p values in the range 
0 < p < n.ka(A/2)? are excited. The bends in one plane do not 
cause coupling to modes with different q values but reduce the delay 
distortion of the modes with different p values. From this physical 
picture we see that the delay distortion problem is reduced to studying 
delay distortion in a slab waveguide. Bends of the fiber in the perpen- 
dicular plane couple the modes with different q but fixed p values. 
Their velocity spread is the same as that considered in the first problem. 
We thus expect to obtain the correct result by considering the delay 
distortion reduction for the mode group with q = 0. 
With g = 0 we obtain from (26) and (57) 


i. (5 7 =) G1 (0)G, (6)d6 


- 8n,A? Uji 2 12(u2 1 ui) 
aCe ee 





Using (24), (36), (54), and (58) we have 
_ 128n2a4A3 0 wu E 12(u? + wu?) i 








e ORQ) | Aa = ue GE = a 
2H4A3 
= 2.26 X 104+. (50) 
The width of the impulse response follows from (52), 
NG — LN 
= OO FRE | et 
at = 0.06" 4 ( i oe) (60) 


From Ref. 9 we obtain the width of the impulse response for uncoupled 
modes 


_ Noli =, 
Ar= De A’, 


(61) 
The improvement that is caused by mode coupling is the ratio of the 
widths of these two impulse responses 


At 0.12a? 


t=; [FQ@La 


(62) 
For the statistical model of a sequence of circularly bent waveguide 
sections we obtain with (48) and (21) 

0.244 ( D )' 


k= ma (63) 





L 
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We may ask for the average radius of curvature that is required to 
cause a ten-fold improvement of the width of the impulse response 
due to mode coupling. Using Z = 1 km and an average length of the 
bent sections of D = 1 cm and the numbers in (50) we obtain for 
R = 0.1, « = 0.022 em— or an average radius of curvature R = 1/x 
= 45 cm. This relatively small radius of curvature may cause very 
substantial excess loss according to the loss example in Section V. 

In order to relate the excess loss to the delay distortion improvement 
we consider the loss penalty that must be paid for a given amount of 
improvement in the width of the impulse response. Both the loss 
formula (40) and the improvement factor R, (62), contain the power 
spectrum of the distortion function f(z). By taking the square of R 
and multiplying it with the loss per length L, o1L, we obtain‘ 


R’o,L = 0.042 = 0.18 dB. (64) 


This important formula is independent of any of the waveguide 
parameters and of the statistics of the axis deformation. This means 
that the loss penalty, o,L, for parabolic index fibers depends only on 
the delay distortion improvement that one wants to achieve. For 
R= 1 we have a loss of oil = 0.18 dB. Clearly, the range of ap- 
plicability of (64) is exceeded in this case since R = 1 means that 
there is no improvement at all. For R = 0.5 we pay a loss penalty of 
0.7 dB, R = 0.1 increases the loss to 18 dB. The already favorable 
delay distortion of the parabolic index fiber can be improved by 
intentional curvature of the waveguide axis. 


VII. DISCUSSION 


We have studied the performance of the parabolic index fiber with 
randomly curved axis. The curvature of the waveguide axis has the 
tendency to force a light beam inside of the fiber towards the fiber 
boundary. In terms of wave optics this means that the wave field 
begins to interact with the boundary of the fiber. If this boundary 
is perfectly smooth no particular harm may be done except that the 
impulse response of the fiber is likely to deteriorate. However, the 
interfaces between two dielectric regions tend to be rough. Surface 
roughness leads to scattering losses. We have thus assumed that the 
interaction of the mode fields with the fiber boundary causes signif- 
icant losses to high-order modes. On this basis we were able to calculate 
the fiber loss caused by random bends of the waveguide axis. For 
bends that approach a sinusoidal shape, with a period comparable to 


t We use o1 of (40) instead of o1 of (42) since R was computed for g = 0. 
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the ray oscillation period in the parabolic index medium, the excess 
losses are extremely high. Bending of the waveguide axis with a period 
equal to the ray oscillation period must be avoided. For a statistical 
model, based on the assumption that the waveguide is composed of a 
sequence of circularly bent sections with random length and random 
radius of curvature, the waveguide losses have been predicted. We 
conclude that average radii of curvature of approximately 1 m can be 
allowed if an excess loss of 10 dB/km can be tolerated. The waveguide 
sections were assumed to have an average length of 1 cm. 

It is possible to reduce the width of the impulse response of a 
parabolic index fiber by coupling its modes by random bends of the 
fiber axis. The impulse response of parabolic index fibers is already 
quite favorable compared to the impulse response of the conventional 
optical fiber with a discontinuous but piecewise constant index dis- 
tribution. Our analysis shows that additional reduction of pulse delay 
distortion is accompanied by losses. A reduction of the pulse width 
to half its uncoupled width increases the loss by 0.7 dB, a ten-fold 
pulse width reduction increases the fiber loss by 18 dB. 
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APPENDIX 


The averaging process used to obtain (42) can be justified as follows. 
Each mode group characterized by the mode number g comprising 


all modes with 
O<p< (./f ncka - a) 


has the loss coefficient o1(q). By definition this can be written 


_ AP) 
oi(q) = P@ 


AP (q) is the power lost from the mode group per unit length and P(q) 
is the power carried by these modes. If we assume, for simplicity, that 
each mode carries the power P we can write P(q) = N(q)P so that we 


have 
AP(q@) 


o1(q) = N(qQ)P ? 





with N(q) indicating the number of modes in the group. The total 
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loss is 


EAP@ = LorlgN(a) 
= ~ P) ~ = N@) 





Replacing the sum by an integral yields formula (42). 
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Effect of Misalignments on Coupling 
Efficiency of Single-Mode Optical 
Fiber Butt Joints 


By J. S. COOK, W. L. MAMMEL, and R. J. GROW 
(Manuscript received April 9, 1973) 


Analysis and computations made here, corroborated by experiment, 
determine the effects of axial displacement and angular misalignment on 
the power coupled between butt-joined, single-mode optical fibers. The 
absolute accuracy with which fibers must be joined on-centers is reduced 
for fibers with relatively smaller core; the angular accuracy 1s increased. 


I. INTRODUCTION 


The lowest-order mode in a clad optical fiber, the hybrid HE 
mode, is the only propagating mode for core sizes less than a few 
wavelengths in diameter (v S$ 2.4).! Hence, small-core, single-mode 
glass fibers are attractive for transmitting optical signals because of 
their potentially low dispersive effects. A possible disadvantage lies 
in the difficulty of joining small-core fibers end-to-end.?* It has been 
suggested‘ that butt joining may be made less critical by reducing the 
size of the fiber in the vicinity of the joint. The computation and ex- 
perimental measurements disclosed here evaluate the advantage to be 
gained by such a procedure. It is found that as the fiber gets smaller, 
the accuracy with which the ends of the fibers to be joined must meet 
on-centers is indeed reduced. At the same time, however, as might be 
expected, the required angular alignment of the fibers becomes more 
critical. 

Both calculations and experiments have been made under the 
assumption that the cladding is of sufficient extent that the role 
played by possible conversion of the zero-order mode to cladding 
modes need not be considered.* The calculation was made by matching 
the fields of the zero-order modes at the joint. 

1439 
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II, ANALYSIS 


The power loss caused by displacement at a joint is readily found 
by first determining the ratio of power accepted by the displaced 
fiber to that presented by the sending fiber, that is, the power-coupling 
ratio. This is determined by assuming that the incident field LE; 
separates at the plane interface (which is perpendicular to the axis 
of the receiving fiber) into the desired zero-order mode that propagates 
in the receiving fiber, and into modes orthogonal to this propagating 
mode. The field of the propagating mode is represented by B-E,, since 
the form £, is known but not the amplitude B. The field of the orthog- 
onal modes is represented by £,. 


i; = B-E, + iy. 


Multiplication by #, and integration over the entire interface gives 
/ E.E,dA = B / E,2dA +0. 
The desired power ratio, c, is then represented by B?. 


gz : E:EpdA / J Baa | 


The power-coupling ratio between fiber ends with parallel but 
translated axes is c:. In order to show the effect of axial displacement 
of a given magnitude, independent of core size, the ratio ci is computed 


TRANSLATION 





Fig. 1—Butt-joined fiber cores displaced by axial translation. 
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as a function of v, the normalized core size,! for several values of the 
parameter d, the normalized displacement. 


y = 78 (m2 — nd), 0) 
d = “(nt — nQy, (2) 


where a is the core radius, s is the axial displacement distance, ) is the 
wavelength, and n, and n, are the core and cladding refractive indexes, 
respectively. 

Figure 1 shows the cross section of a fiber displaced by axial trans- 
lation. R and S are normalized to radius a. The origin is defined as 
the midpoint between core centers. R and @ are the polar coordinates 
of an arbitrary point P in the interface. R, is the distance from P to 
the center of the first fiber, and Rz is the distance to the center of the 
second fiber. 

S 


2 
R? = Re+ (5) — RS cos 8, 


S 2 
R3= R+ (5) + RS cos 0. 


Let Ai be the area of the interface of the displaced fibers, and A 
be the area of the cross section of the sending fiber ; then 


a = | [ _B(R)E(R:) dA f [ (RA | 


The function E£ ds defined! by 


2 





J o(uR) 
Tay? fi SL, 
ee oe oR) 
o\W. 
“K,(w) , R> 1, 


where J, and K, are the regular and modified Bessel functions of 
order zero. 
The values of u and w are determined from the eigenvalue equations, 


v= (w+ w’)), 
ud 1(u) = wKi(w) 
Jo(u)  - K(w) 


The integral in the denominator of the expression ¢c: can be inte- 
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grated analytically. For an infinite cladding, 


feos = of S75 | 





For a fiber of radius R,, (R. > 1), 


a E(R)dA = (2a) = a [K3(wR.) — K2(wR.)]. 


The integral in the numerator can be divided so that portions of the 
integration can be done analytically. 

The power-coupling ratio between fibers with angular displacement 
of the fiber axes is ce. In order to show the effect of a fixed angular 
displacement of the axes for different core sizes, the ratio cz is com- 
puted as a function of v for several values of the parameter b, the 
normalized displacement angle, 


b= sin ¢ Be. 
V1 — n2/n%— V2A’ 


where ¢ is the angle between the fiber axes, and A is the ratio of core- 
to-cladding index difference to the core index. 

Figure 2 shows a cross section of the fiber joint in the plane of the 
axes of the fiber and auxiliary cross sections perpendicular to the axes 
of the fibers. R’, 6’, 2’ and 2’, y’, 2’ are coordinate systems oriented 
with respect to the sending fiber while R, 0, z and 2, y, 2 are coordinate 
systems oriented with respect to the receiving fiber; ¢@ is the angle 
between the axes of the two fibers (and therefore between the planes 
perpendicular to the axes). From Fig. 2 it can be seen that 


R’ = R(1 — sin’6 sin’¢)}. 
The field of the sending fiber at the point P of the interface Ly is 
E;(R) = E(B’) cos 62’, 


where @ is the normalized propagation constant. 


(3) 





v 

ae 
For the small angles under consideration the maximum deviation of 
R’/R from 1 is at most 2-10-?, so it is feasible to approximate R’ by R. 
Therefore, 62’ = bvR sin 6 and 
2 


eo(0) = J, E*(R) cos (bvR sin 6)dA i ie (Ryd A | 
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CLADDING 


PIN TERGAGE 





Fig. 2—Analytical model of butt-joined fiber cores displaced by angular mis- 
alignment, ¢. 


The integral in the numerator, since we have approximated R’ by R, 
can be converted to a single integral. 


Re Qr Re 
i i, E?(R) cos (buR sin 6)RdedR = 2n i E?(R)J9(bvR) RAR. 
0 0 


The integral in the denominator is the same as before and can be 
integrated analytically. 


III. EXPERIMENT 


Experimental verification was carried out at microwave frequencies 
because equipment was readily available and dimensional control 
more certain than at optical frequencies. Polyfoam served as the core, 
and air as the cladding. The experimental arrangement is shown in 
Fig. 8. The polyfoam rod had dielectric constant of 1.06; hence, 
(n2 — n?)* = 0.245. Results of computations from the analysis and 
experimentally measured points are shown overlaid in Figs. 4 and 5. 
The vertical point dimension indicates the disparity between two 
systematic measurements. Disparity between the analytical and 
experimental results for very small core dimensions and larger coupling 
loss is not understood; but it is not critical to the conclusions. 

Figure 6 shows an overlay of calculated power coupling as a function 
of v due to both angular and lateral misalignments. It can be argued 
that the total power lost at the junction is roughly equal to the sum 
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Fig. 3—Experimental arrangement used to corroborate analysis using microwave-guiding polyfoam rods. (a) Coupling vs axial 
displacement measurement configuration. (b) Coupling vs angular misalignment measurement configuration. 
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Fig. 4—Power coupling through translationally displaced joint as a function of 
normalized core size and axial translation, d. 
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Fig. 5—Power coupling through angularly misaligned joint as a function of normal- 
ized core size and misalignment, b. 
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Fig. 6—Minimum power coupling through example joint as a function of normalized 
core size. Dashed line shows power coupled at maximum axial translation, dotted 
line shows composite minimum coupling. 


of the lost power due to the two kinds of misalignment. Since many 
modes are involved in the scattered energy, and the displacement 
coordinates are different, the scattered modes must be essentially 
orthogonal, hence, power-additive. The dash-dot curves in Fig. 6 
show the fraction of power lying within 10 and 20 radii of the core 
as labeled. 


IV. CONCLUSIONS 


In general, considering the nature of field spreading that accom- 
panies the decrease in single-mode fiber core diameter, one concludes 
what one would expect. As the core size parameter, v, decreases below 
about 2, the fields spread, and axial displacement at the joints becomes 
less and less critical since more and more overlap of fields results from 
a given offset. At the same time, the effective “aperture” at the fiber 
end increases, the “antenna”’ becomes more directive, and the angular 
alignment becomes more critical. 

More detailed conclusions must be drawn from the particular 
joining problem at hand. Suppose, for example, one wants to join 
single-mode fibers of 3-mil diameter, having a core loaded to produce 
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an index difference at the interface of 4 percent. If the refractive index 
of the fiber material is about 1.5, then 


(n2 — n2)t ey n.V2A = 0.15. 


If the fiber is designed to have v = 2.0 and the optical wavelength, \, 
is 1 um, the core radius may be found from eq. (1) to be 


ary 2.1 pm. 


Since 1.5 mils & 38 um, the cladding radius is about 18 times the core 
radius, so 


Rk. = 18. 


Suppose the centering accuracy of the tool or fixture that will be 
used in aligning the fiber ends at a joint is about +1 um. (It pre- 
sumably will center the fiber with reference to its o.d.) And suppose 
the core is centered in the fiber with +1-percent accuracy, that is, 
the core center is never displaced from the fiber center by more than 
1 percent of the fiber diameter. The net maximum displacement, 
then, will be 

Sn A 175 ph 
and, from (2), 
dm & 1.65. 


If the fiber is drawn down to a smaller size at the end to decrease 
the lateral displacement sensitivity at the joint, it is reasonable to 
assume that the alignment tool maintains the same l-um accuracy 
and the core maintains the same +1-percent centering accuracy; the 
latter of which produces a net decrease in dm. The dashed curve in 
Fig. 6 shows what happens to the minimum coupling, ci, as v scales 
down with fiber size. 

If at the same time one assumes that the net angular misalignment, 
¢, is less than about 0.01 radian, then from (3), 


b 0.1, 


and the combined effect of angular and lateral displacement varies 
with v as shown by the dotted line in Fig. 6. 

The total worst-case joint loss in this example, then, may be reduced 
by about a factor of two by drawing the fiber (core and cladding 
together) down to half its normal size at the ends for joining. At 
v = 1 the field extending outside the cladding is still negligible. 
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